execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

ex5 regression
ex5.stan

data{
  int N;
  int K;
  vector[N] y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
  real<lower=0> s;
}
model{
  vector[N] m=X*b;
  y~normal(m,s);
}
generated quantities{
  vector[N] y1;
  vector[N] m1=X*b;
  for(i in 1:N){
    y1[i]=normal_rng(m1[i],s);
  }
}

normal linear models

n=30
mdl=cmdstan_model('./ex5.stan')

single regression

tb=tibble(x=runif(n,0,9),
          y=rnorm(n,x,1))
f=formula(y~x)
par(mfrow=c(1,1))
plot(tb$x,tb$y)

qplot(data=tb,x,y)

X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -105.796 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       15      -12.8956   0.000206144   3.95095e-05           1           1       19    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
mle
##  variable estimate
##    lp__     -12.90
##    b[1]      -0.20
##    b[2]       0.99
##    s          0.93
##    y1[1]     -0.83
##    y1[2]      8.06
##    y1[3]      1.77
##    y1[4]     -0.74
##    y1[5]      7.03
##    y1[6]      7.43
##    y1[7]      4.32
##    y1[8]      6.94
##    y1[9]      6.45
##    y1[10]     3.48
##    y1[11]     8.73
##    y1[12]     6.94
##    y1[13]     0.73
##    y1[14]     7.45
##    y1[15]     5.16
##    y1[16]     3.37
##    y1[17]    -0.48
##    y1[18]     0.65
##    y1[19]     0.76
##    y1[20]     6.17
##    y1[21]     6.42
##    y1[22]     1.03
##    y1[23]     8.17
##    y1[24]     5.51
##    y1[25]     8.03
##    y1[26]     3.09
##    y1[27]     7.77
##    y1[28]     0.37
##    y1[29]     8.27
##    y1[30]     2.45
##    m1[1]      0.39
##    m1[2]      6.91
##    m1[3]      3.32
##    m1[4]      0.28
##    m1[5]      6.06
##    m1[6]      7.26
##    m1[7]      6.42
##    m1[8]      5.90
##    m1[9]      5.79
##    m1[10]     4.45
##    m1[11]     7.98
##    m1[12]     6.72
##    m1[13]     1.55
##    m1[14]     7.37
##    m1[15]     6.36
##    m1[16]     3.20
##    m1[17]     0.00
##    m1[18]     0.81
##    m1[19]     2.79
##    m1[20]     5.89
##    m1[21]     7.43
##    m1[22]     1.11
##    m1[23]     8.03
##    m1[24]     6.23
##    m1[25]     6.85
##    m1[26]     1.86
##    m1[27]     7.50
##    m1[28]     1.12
##    m1[29]     8.31
##    m1[30]     1.70
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -14.54 -14.21 1.26 1.05 -17.07 -13.13 1.00      718      868
##    b[1]    -0.23  -0.23 0.37 0.38  -0.83   0.35 1.01      391      386
##    b[2]     1.00   1.00 0.07 0.07   0.89   1.11 1.00      539      749
##    s        1.01   0.99 0.14 0.13   0.80   1.26 1.00     1110      772
##    y1[1]    0.38   0.40 1.05 1.01  -1.30   2.11 1.00     1188     1852
##    y1[2]    6.88   6.88 1.07 1.03   5.13   8.66 1.00     1677     1742
##    y1[3]    3.33   3.31 1.02 1.01   1.64   5.02 1.00     1641     1745
##    y1[4]    0.21   0.20 1.07 1.05  -1.58   1.96 1.00     1662     1809
##    y1[5]    6.04   6.04 1.02 1.00   4.31   7.70 1.00     2009     1801
##    y1[6]    7.26   7.23 1.03 1.04   5.67   8.91 1.00     1994     2058
##    y1[7]    6.42   6.44 1.06 1.04   4.60   8.13 1.00     1718     1860
##    y1[8]    5.92   5.92 1.03 1.01   4.23   7.59 1.00     2131     2009
##    y1[9]    5.78   5.78 0.99 0.94   4.19   7.43 1.00     2115     2015
##    y1[10]   4.44   4.45 1.07 1.07   2.67   6.14 1.00     2023     1972
##    y1[11]   7.99   7.98 1.03 0.98   6.28   9.69 1.00     2061     1968
##    y1[12]   6.71   6.71 1.05 1.02   4.98   8.39 1.00     1853     1853
##    y1[13]   1.52   1.55 1.06 1.05  -0.26   3.23 1.00     1906     1932
##    y1[14]   7.39   7.36 1.09 1.05   5.59   9.18 1.00     2039     1684
##    y1[15]   6.36   6.37 1.03 0.98   4.64   8.04 1.00     2102     1927
##    y1[16]   3.18   3.16 1.05 1.03   1.45   4.91 1.00     1814     1734
##    y1[17]   0.03   0.05 1.09 1.09  -1.77   1.83 1.00     1573     1888
##    y1[18]   0.82   0.88 1.12 1.08  -1.08   2.61 1.00     1208     1622
##    y1[19]   2.78   2.76 1.04 0.99   1.14   4.55 1.00     2016     1950
##    y1[20]   5.90   5.89 1.05 1.01   4.18   7.64 1.00     1991     1856
##    y1[21]   7.47   7.45 1.09 1.05   5.67   9.21 1.00     1949     1881
##    y1[22]   1.12   1.12 1.07 1.05  -0.58   2.89 1.00     1722     1820
##    y1[23]   8.06   8.05 1.07 1.02   6.28   9.78 1.00     2072     1758
##    y1[24]   6.24   6.26 1.04 1.02   4.53   7.90 1.00     1963     1997
##    y1[25]   6.82   6.82 1.05 1.06   5.11   8.57 1.00     1987     1748
##    y1[26]   1.80   1.80 1.03 1.00   0.12   3.44 1.00     1819     1827
##    y1[27]   7.49   7.51 1.06 1.04   5.71   9.20 1.00     2016     2010
##    y1[28]   1.10   1.07 1.07 1.09  -0.57   2.89 1.00     1583     1806
##    y1[29]   8.32   8.31 1.07 1.06   6.63  10.09 1.00     2110     2014
##    y1[30]   1.71   1.70 1.03 1.02   0.07   3.39 1.00     1738     1968
##    m1[1]    0.37   0.37 0.34 0.34  -0.18   0.91 1.01      389      394
##    m1[2]    6.92   6.92 0.24 0.24   6.52   7.32 1.00     2296     1605
##    m1[3]    3.31   3.31 0.20 0.20   2.98   3.65 1.00      508      985
##    m1[4]    0.25   0.26 0.34 0.35  -0.30   0.80 1.01      389      394
##    m1[5]    6.06   6.06 0.21 0.21   5.71   6.41 1.00     2643     1606
##    m1[6]    7.27   7.26 0.26 0.25   6.84   7.70 1.00     2083     1647
##    m1[7]    6.42   6.42 0.22 0.22   6.05   6.79 1.00     2611     1581
##    m1[8]    5.91   5.91 0.21 0.20   5.57   6.25 1.00     2566     1606
##    m1[9]    5.79   5.79 0.20 0.20   5.46   6.13 1.00     2478     1665
##    m1[10]   4.44   4.44 0.19 0.18   4.13   4.74 1.00      878     1480
##    m1[11]   8.00   7.99 0.30 0.29   7.51   8.49 1.00     1670     1591
##    m1[12]   6.72   6.72 0.24 0.23   6.33   7.11 1.00     2433     1616
##    m1[13]   1.53   1.53 0.28 0.28   1.08   1.97 1.01      396      641
##    m1[14]   7.39   7.38 0.27 0.26   6.95   7.83 1.00     2006     1619
##    m1[15]   6.36   6.36 0.22 0.22   5.99   6.73 1.00     2634     1581
##    m1[16]   3.19   3.19 0.21 0.21   2.85   3.53 1.00      491      932
##    m1[17]  -0.03  -0.02 0.36 0.36  -0.61   0.54 1.01      390      386
##    m1[18]   0.78   0.78 0.31 0.32   0.27   1.29 1.01      389      593
##    m1[19]   2.77   2.77 0.22 0.22   2.41   3.13 1.01      449      773
##    m1[20]   5.90   5.90 0.21 0.20   5.56   6.24 1.00     2562     1609
##    m1[21]   7.44   7.43 0.27 0.26   7.00   7.88 1.00     1973     1619
##    m1[22]   1.09   1.09 0.30 0.30   0.61   1.57 1.01      390      604
##    m1[23]   8.04   8.04 0.30 0.29   7.55   8.54 1.00     1651     1554
##    m1[24]   6.24   6.23 0.22 0.21   5.87   6.59 1.00     2665     1607
##    m1[25]   6.86   6.85 0.24 0.23   6.45   7.25 1.00     2333     1616
##    m1[26]   1.84   1.84 0.26 0.26   1.42   2.26 1.01      402      683
##    m1[27]   7.51   7.51 0.27 0.26   7.07   7.96 1.00     1929     1603
##    m1[28]   1.10   1.10 0.30 0.30   0.62   1.58 1.01      390      604
##    m1[29]   8.32   8.32 0.31 0.30   7.80   8.84 1.00     1531     1487
##    m1[30]   1.68   1.68 0.27 0.27   1.24   2.11 1.01      398      673

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)

par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(1:n,y,data=tb,col=I('red'))+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

multiple regression

tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
          y=rnorm(n,x1-x2,1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -15999 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       44      -13.1144   0.000133986   0.000327327           1           1       86    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -13.11
##    b[1]      -0.16
##    b[2]       1.04
##    b[3]      -0.95
##    s          0.94
##    y1[1]      0.64
##    y1[2]      2.63
##    y1[3]     -3.29
##    y1[4]      2.85
##    y1[5]     -3.74
##    y1[6]      4.86
##    y1[7]     -2.92
##    y1[8]     -2.71
##    y1[9]      8.13
##    y1[10]     1.65
##    y1[11]    -1.58
##    y1[12]     1.49
##    y1[13]    -4.53
##    y1[14]    -0.26
##    y1[15]    -1.97
##    y1[16]     2.93
##    y1[17]    -0.68
##    y1[18]     2.83
##    y1[19]     1.95
##    y1[20]    -1.32
##    y1[21]    -0.38
##    y1[22]    -5.81
##    y1[23]    -1.32
##    y1[24]    -2.67
##    y1[25]     7.60
##    y1[26]     3.43
##    y1[27]     2.96
##    y1[28]     5.53
##    y1[29]    -4.26
##    y1[30]     5.03
##    m1[1]      1.10
##    m1[2]      1.35
##    m1[3]     -3.59
##    m1[4]      3.95
##    m1[5]     -3.21
##    m1[6]      4.05
##    m1[7]     -2.41
##    m1[8]     -3.91
##    m1[9]      7.81
##    m1[10]     1.16
##    m1[11]    -0.61
##    m1[12]     3.02
##    m1[13]    -4.30
##    m1[14]    -0.16
##    m1[15]    -3.05
##    m1[16]     2.30
##    m1[17]    -0.76
##    m1[18]     2.14
##    m1[19]     2.73
##    m1[20]    -2.01
##    m1[21]    -1.38
##    m1[22]    -4.31
##    m1[23]    -1.23
##    m1[24]    -4.30
##    m1[25]     5.22
##    m1[26]     3.06
##    m1[27]     2.51
##    m1[28]     4.29
##    m1[29]    -3.01
##    m1[30]     4.10
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -15.38 -15.04 1.52 1.32 -18.45 -13.61 1.00      623      644
##    b[1]    -0.10  -0.13 0.54 0.53  -0.98   0.84 1.01      330      317
##    b[2]     1.03   1.04 0.09 0.09   0.88   1.18 1.00      547      555
##    b[3]    -0.96  -0.96 0.08 0.07  -1.08  -0.83 1.01      674      710
##    s        1.04   1.02 0.16 0.15   0.82   1.32 1.00     1072      774
##    y1[1]    1.09   1.08 1.13 1.12  -0.78   2.93 1.00     1109     1268
##    y1[2]    1.33   1.30 1.11 1.07  -0.55   3.18 1.00     1984     1784
##    y1[3]   -3.52  -3.52 1.09 1.06  -5.31  -1.76 1.00     1941     1835
##    y1[4]    3.91   3.91 1.12 1.10   2.07   5.74 1.00     1615     2038
##    y1[5]   -3.19  -3.18 1.08 1.03  -4.98  -1.40 1.00     1928     1668
##    y1[6]    4.02   4.05 1.09 1.06   2.20   5.73 1.00     1728     1656
##    y1[7]   -2.35  -2.36 1.09 1.05  -4.13  -0.49 1.00     1677     1548
##    y1[8]   -3.91  -3.93 1.13 1.10  -5.74  -2.01 1.00     1975     1994
##    y1[9]    7.77   7.77 1.18 1.12   5.87   9.73 1.00     2031     1853
##    y1[10]   1.16   1.17 1.07 1.05  -0.64   2.89 1.00     1885     1961
##    y1[11]  -0.62  -0.59 1.08 1.01  -2.41   1.13 1.00     1878     1838
##    y1[12]   3.06   3.04 1.13 1.10   1.22   4.91 1.00     1869     1889
##    y1[13]  -4.24  -4.27 1.11 1.05  -6.03  -2.40 1.00     2000     2015
##    y1[14]  -0.12  -0.10 1.07 1.04  -1.84   1.65 1.00     1769     1822
##    y1[15]  -3.04  -3.05 1.12 1.12  -4.85  -1.27 1.00     1748     1754
##    y1[16]   2.27   2.26 1.14 1.09   0.40   4.10 1.00     1818     1637
##    y1[17]  -0.77  -0.75 1.07 1.09  -2.51   0.94 1.00     1934     1487
##    y1[18]   2.15   2.12 1.05 1.04   0.42   3.86 1.00     1940     1930
##    y1[19]   2.68   2.66 1.09 1.07   0.95   4.49 1.00     1871     1665
##    y1[20]  -2.00  -1.96 1.11 1.08  -3.85  -0.23 1.00     1765     1767
##    y1[21]  -1.38  -1.40 1.07 1.06  -3.15   0.34 1.00     1925     1770
##    y1[22]  -4.30  -4.30 1.11 1.11  -6.07  -2.43 1.00     2100     1932
##    y1[23]  -1.13  -1.18 1.17 1.15  -3.04   0.84 1.00     1414     1401
##    y1[24]  -4.30  -4.33 1.11 1.09  -6.05  -2.40 1.00     1508     1915
##    y1[25]   5.23   5.21 1.13 1.08   3.37   7.10 1.00     1760     1811
##    y1[26]   3.01   3.02 1.09 1.05   1.22   4.75 1.00     1957     1791
##    y1[27]   2.52   2.48 1.10 1.09   0.77   4.39 1.00     1934     1830
##    y1[28]   4.31   4.29 1.14 1.13   2.46   6.16 1.00     1909     1760
##    y1[29]  -3.01  -3.01 1.10 1.05  -4.84  -1.25 1.00     2071     1973
##    y1[30]   4.12   4.12 1.14 1.08   2.26   5.94 1.00     2159     1921
##    m1[1]    1.13   1.13 0.42 0.40   0.45   1.85 1.01      343      360
##    m1[2]    1.37   1.36 0.26 0.25   0.96   1.79 1.00      485      610
##    m1[3]   -3.56  -3.56 0.31 0.29  -4.06  -3.05 1.00      947     1374
##    m1[4]    3.92   3.93 0.44 0.43   3.16   4.62 1.00      734      836
##    m1[5]   -3.19  -3.19 0.29 0.27  -3.65  -2.71 1.00      957     1414
##    m1[6]    4.05   4.05 0.30 0.29   3.57   4.56 1.00     1766     1382
##    m1[7]   -2.37  -2.38 0.35 0.34  -2.94  -1.79 1.01      427      454
##    m1[8]   -3.91  -3.91 0.37 0.37  -4.51  -3.30 1.00     1240     1227
##    m1[9]    7.80   7.78 0.48 0.47   7.01   8.60 1.00     1818     1254
##    m1[10]   1.16   1.16 0.20 0.20   0.84   1.50 1.00     1963     1537
##    m1[11]  -0.61  -0.61 0.29 0.30  -1.11  -0.14 1.00      948     1303
##    m1[12]   3.03   3.03 0.25 0.24   2.62   3.46 1.00     1571     1248
##    m1[13]  -4.27  -4.27 0.37 0.36  -4.87  -3.64 1.00      695      791
##    m1[14]  -0.17  -0.17 0.27 0.27  -0.64   0.26 1.00     1069     1340
##    m1[15]  -3.06  -3.06 0.41 0.40  -3.72  -2.40 1.00      821     1026
##    m1[16]   2.27   2.28 0.38 0.37   1.62   2.86 1.00      711      787
##    m1[17]  -0.76  -0.76 0.22 0.22  -1.11  -0.39 1.00     2137     1576
##    m1[18]   2.15   2.14 0.23 0.23   1.79   2.54 1.00      947     1086
##    m1[19]   2.71   2.71 0.33 0.33   2.16   3.26 1.00      988     1100
##    m1[20]  -1.97  -1.98 0.39 0.38  -2.60  -1.31 1.01      379      356
##    m1[21]  -1.36  -1.36 0.24 0.24  -1.74  -0.96 1.00      615      835
##    m1[22]  -4.29  -4.30 0.33 0.33  -4.83  -3.75 1.00     1755     1611
##    m1[23]  -1.19  -1.20 0.46 0.45  -1.93  -0.41 1.01      346      332
##    m1[24]  -4.29  -4.30 0.36 0.37  -4.87  -3.69 1.00     1511     1290
##    m1[25]   5.22   5.21 0.37 0.36   4.64   5.83 1.00     1027      996
##    m1[26]   3.04   3.04 0.30 0.30   2.55   3.55 1.00     1579     1560
##    m1[27]   2.52   2.52 0.29 0.28   2.05   3.00 1.00      522      725
##    m1[28]   4.29   4.29 0.32 0.31   3.79   4.82 1.00     1123     1085
##    m1[29]  -3.01  -3.02 0.34 0.34  -3.55  -2.45 1.00     1130     1130
##    m1[30]   4.10   4.09 0.30 0.29   3.61   4.61 1.00     1181     1205

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(1:n,y,data=tb,col=I('red'))+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

ANOVA

tb=tibble(c=sample(c('a','b','c'),n,replace=T),
          y=rnorm(n,(c=='b')*2-(c=='c')*2,1))
f=formula(y~c)
par(mfrow=c(1,1))
boxplot(y~c,tb)

qplot(data=tb,c,y,geom='boxplot')

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -53.587 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       17      -9.48161    0.00018107   0.000101393           1           1       26    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__      -9.48
##    b[1]       0.41
##    b[2]       1.02
##    b[3]      -2.56
##    s          0.83
##    y1[1]     -0.72
##    y1[2]     -0.37
##    y1[3]      1.08
##    y1[4]     -3.57
##    y1[5]     -2.43
##    y1[6]      0.66
##    y1[7]      0.70
##    y1[8]     -2.74
##    y1[9]      1.29
##    y1[10]     0.22
##    y1[11]     0.52
##    y1[12]    -1.69
##    y1[13]     0.10
##    y1[14]    -3.04
##    y1[15]     0.34
##    y1[16]     0.81
##    y1[17]     0.63
##    y1[18]    -0.22
##    y1[19]    -0.24
##    y1[20]     0.78
##    y1[21]     0.59
##    y1[22]     1.58
##    y1[23]    -0.75
##    y1[24]    -0.19
##    y1[25]    -1.82
##    y1[26]     0.39
##    y1[27]     0.53
##    y1[28]     1.48
##    y1[29]    -0.08
##    y1[30]     0.84
##    m1[1]      0.41
##    m1[2]      0.41
##    m1[3]      0.41
##    m1[4]     -2.15
##    m1[5]     -2.15
##    m1[6]      1.43
##    m1[7]      1.43
##    m1[8]     -2.15
##    m1[9]      1.43
##    m1[10]     0.41
##    m1[11]     0.41
##    m1[12]     0.41
##    m1[13]     0.41
##    m1[14]    -2.15
##    m1[15]     1.43
##    m1[16]     1.43
##    m1[17]     1.43
##    m1[18]     0.41
##    m1[19]     0.41
##    m1[20]     1.43
##    m1[21]     0.41
##    m1[22]     0.41
##    m1[23]     0.41
##    m1[24]    -2.15
##    m1[25]    -2.15
##    m1[26]     0.41
##    m1[27]     0.41
##    m1[28]     1.43
##    m1[29]     0.41
##    m1[30]     1.43
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -11.80 -11.50 1.49 1.27 -14.65 -10.07 1.00      829      937
##    b[1]     0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    b[2]     1.00   1.00 0.39 0.40   0.39   1.65 1.00      901     1041
##    b[3]    -2.54  -2.55 0.46 0.46  -3.26  -1.78 1.01      780      994
##    s        0.92   0.90 0.13 0.13   0.73   1.16 1.00     1829     1485
##    y1[1]    0.40   0.38 0.98 0.95  -1.20   1.98 1.00     1819     1893
##    y1[2]    0.41   0.43 0.98 0.97  -1.18   1.97 1.00     2033     1810
##    y1[3]    0.44   0.44 0.97 0.94  -1.18   2.03 1.00     1795     2090
##    y1[4]   -2.10  -2.10 0.99 0.97  -3.68  -0.43 1.00     2123     1929
##    y1[5]   -2.14  -2.14 1.01 1.00  -3.76  -0.49 1.00     1806     1853
##    y1[6]    1.42   1.38 0.97 0.95  -0.14   3.10 1.00     1820     1722
##    y1[7]    1.42   1.42 0.97 0.94  -0.23   3.02 1.01     1705     1632
##    y1[8]   -2.07  -2.08 1.01 0.96  -3.71  -0.40 1.00     1860     1971
##    y1[9]    1.39   1.36 0.97 0.94  -0.21   2.99 1.00     1878     2015
##    y1[10]   0.42   0.44 0.93 0.90  -1.14   1.92 1.00     1898     1971
##    y1[11]   0.44   0.42 0.94 0.92  -1.10   1.99 1.00     1948     2001
##    y1[12]   0.38   0.40 0.95 0.96  -1.16   1.90 1.00     1553     1959
##    y1[13]   0.42   0.40 0.97 0.95  -1.12   2.04 1.00     1992     1943
##    y1[14]  -2.13  -2.13 0.98 0.96  -3.68  -0.54 1.00     1716     1868
##    y1[15]   1.41   1.41 0.98 0.96  -0.17   2.99 1.01     1984     1824
##    y1[16]   1.45   1.48 0.96 0.88  -0.18   3.01 1.00     2149     1973
##    y1[17]   1.40   1.42 0.99 0.98  -0.17   3.04 1.00     1961     2015
##    y1[18]   0.42   0.41 0.95 0.93  -1.14   1.99 1.00     1872     1792
##    y1[19]   0.43   0.40 0.95 0.91  -1.09   2.03 1.00     1885     1748
##    y1[20]   1.40   1.40 0.98 0.90  -0.22   3.01 1.00     1920     1885
##    y1[21]   0.40   0.42 0.97 0.94  -1.20   2.00 1.00     1710     1784
##    y1[22]   0.42   0.40 0.96 0.92  -1.17   2.00 1.00     1980     1891
##    y1[23]   0.38   0.38 0.98 0.98  -1.24   1.98 1.00     1892     1740
##    y1[24]  -2.11  -2.14 1.01 0.98  -3.79  -0.43 1.00     1887     1991
##    y1[25]  -2.15  -2.16 1.01 1.01  -3.78  -0.44 1.00     1672     1842
##    y1[26]   0.42   0.43 0.94 0.95  -1.12   1.95 1.00     1821     1757
##    y1[27]   0.41   0.40 0.97 0.93  -1.14   2.00 1.00     1807     1915
##    y1[28]   1.41   1.39 0.98 0.96  -0.20   2.95 1.00     1575     1934
##    y1[29]   0.38   0.37 0.96 0.94  -1.17   1.94 1.00     1959     1912
##    y1[30]   1.40   1.41 0.99 0.91  -0.24   3.00 1.00     1917     1869
##    m1[1]    0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[2]    0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[3]    0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[4]   -2.13  -2.15 0.38 0.37  -2.74  -1.50 1.00     1103     1304
##    m1[5]   -2.13  -2.15 0.38 0.37  -2.74  -1.50 1.00     1103     1304
##    m1[6]    1.42   1.42 0.31 0.31   0.92   1.91 1.00     1442     1591
##    m1[7]    1.42   1.42 0.31 0.31   0.92   1.91 1.00     1442     1591
##    m1[8]   -2.13  -2.15 0.38 0.37  -2.74  -1.50 1.00     1103     1304
##    m1[9]    1.42   1.42 0.31 0.31   0.92   1.91 1.00     1442     1591
##    m1[10]   0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[11]   0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[12]   0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[13]   0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[14]  -2.13  -2.15 0.38 0.37  -2.74  -1.50 1.00     1103     1304
##    m1[15]   1.42   1.42 0.31 0.31   0.92   1.91 1.00     1442     1591
##    m1[16]   1.42   1.42 0.31 0.31   0.92   1.91 1.00     1442     1591
##    m1[17]   1.42   1.42 0.31 0.31   0.92   1.91 1.00     1442     1591
##    m1[18]   0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[19]   0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[20]   1.42   1.42 0.31 0.31   0.92   1.91 1.00     1442     1591
##    m1[21]   0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[22]   0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[23]   0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[24]  -2.13  -2.15 0.38 0.37  -2.74  -1.50 1.00     1103     1304
##    m1[25]  -2.13  -2.15 0.38 0.37  -2.74  -1.50 1.00     1103     1304
##    m1[26]   0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[27]   0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[28]   1.42   1.42 0.31 0.31   0.92   1.91 1.00     1442     1591
##    m1[29]   0.41   0.42 0.24 0.23   0.01   0.80 1.00      924     1376
##    m1[30]   1.42   1.42 0.31 0.31   0.92   1.91 1.00     1442     1591

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')

qplot(data=tb,y,y1,shape=c,size=I(2),xlab='obs.',ylab='prd.')

plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')

qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(data=tb,1:n,y,col=c)+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

ANCOVA

tb=tibble(x=runif(n,0,9),c=sample(c('a','b','c'),n,replace=T),
          y=rnorm(n,x+(c=='b')*2-(c=='c')*2,1))

f=formula(y~x+c)
par(mfrow=c(1,1))
#plot(tb$x1,tb$y,pch=tb$c1)

lv=c(0,1,2)
plot(tb$x,tb$y,pch=lv[factor(tb$c)])

qplot(data=tb,x,y,shape=c,size=I(2))

plot(tb$x,tb$y,col=1+lv[factor(tb$c)])

qplot(data=tb,x,y,col=c)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
## Initial log joint probability = -264.755 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       42      -10.8678   6.35115e-05   0.000874227      0.8883      0.8883       58    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -10.87
##    b[1]       0.41
##    b[2]       0.97
##    b[3]       1.95
##    b[4]      -1.61
##    s          0.87
##    y1[1]      7.29
##    y1[2]      5.29
##    y1[3]      5.51
##    y1[4]      7.13
##    y1[5]      4.15
##    y1[6]      6.09
##    y1[7]     -0.12
##    y1[8]      5.09
##    y1[9]      5.60
##    y1[10]     2.20
##    y1[11]     9.16
##    y1[12]     7.63
##    y1[13]     3.24
##    y1[14]     9.70
##    y1[15]     4.27
##    y1[16]     5.05
##    y1[17]     0.78
##    y1[18]     7.85
##    y1[19]     6.65
##    y1[20]    10.14
##    y1[21]     8.69
##    y1[22]     0.07
##    y1[23]     8.23
##    y1[24]     3.84
##    y1[25]     4.40
##    y1[26]     9.11
##    y1[27]    10.32
##    y1[28]     3.01
##    y1[29]    11.96
##    y1[30]     4.94
##    m1[1]      5.16
##    m1[2]      4.68
##    m1[3]      7.18
##    m1[4]      6.21
##    m1[5]      3.16
##    m1[6]      5.42
##    m1[7]      0.90
##    m1[8]      5.45
##    m1[9]      6.15
##    m1[10]     2.63
##    m1[11]     8.86
##    m1[12]     8.13
##    m1[13]     3.83
##    m1[14]     9.75
##    m1[15]     3.94
##    m1[16]     5.96
##    m1[17]     0.55
##    m1[18]     7.32
##    m1[19]     7.38
##    m1[20]     9.94
##    m1[21]     9.64
##    m1[22]     0.26
##    m1[23]     7.39
##    m1[24]     3.09
##    m1[25]     4.46
##    m1[26]     9.23
##    m1[27]    10.20
##    m1[28]     1.08
##    m1[29]    10.42
##    m1[30]     5.88
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='m',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -13.71 -13.37 1.70 1.50 -16.92 -11.58 1.00      752     1036
##    b[1]     0.45   0.46 0.41 0.40  -0.23   1.11 1.00      863      945
##    b[2]     0.97   0.97 0.08 0.08   0.84   1.09 1.00     1275     1273
##    b[3]     1.95   1.94 0.43 0.43   1.24   2.66 1.00      854      929
##    b[4]    -1.60  -1.64 0.58 0.57  -2.50  -0.61 1.00      640      836
##    s        0.98   0.97 0.14 0.13   0.78   1.24 1.00     1526     1177
##    y1[1]    5.16   5.14 1.10 1.09   3.31   6.98 1.00     1901     1868
##    y1[2]    4.67   4.68 1.03 1.00   2.99   6.31 1.00     1739     1895
##    y1[3]    7.13   7.13 1.07 1.03   5.38   8.93 1.00     1814     1936
##    y1[4]    6.17   6.19 1.12 1.07   4.32   7.98 1.00     1788     1848
##    y1[5]    3.14   3.13 1.02 0.97   1.51   4.83 1.00     1685     1830
##    y1[6]    5.39   5.41 1.06 1.03   3.62   7.07 1.00     1704     1740
##    y1[7]    0.93   0.94 1.05 0.99  -0.82   2.64 1.00     1455     1724
##    y1[8]    5.45   5.45 1.08 1.06   3.73   7.25 1.00     1698     1820
##    y1[9]    6.12   6.09 1.02 0.96   4.46   7.82 1.00     1946     1853
##    y1[10]   2.69   2.71 1.03 1.03   0.98   4.34 1.00     1851     1676
##    y1[11]   8.86   8.82 1.03 1.00   7.17  10.61 1.00     1816     2012
##    y1[12]   8.15   8.13 0.98 0.96   6.53   9.77 1.00     2015     1981
##    y1[13]   3.86   3.85 1.09 1.03   2.15   5.70 1.00     1831     1859
##    y1[14]   9.72   9.75 1.03 0.97   8.04  11.40 1.00     2052     1915
##    y1[15]   3.97   3.95 1.04 1.02   2.29   5.66 1.00     1807     1889
##    y1[16]   5.93   5.91 1.00 1.01   4.24   7.54 1.00     1941     1801
##    y1[17]   0.55   0.53 1.04 0.97  -1.13   2.31 1.00     1378     1802
##    y1[18]   7.32   7.30 1.00 1.02   5.68   8.89 1.00     2082     1743
##    y1[19]   7.38   7.38 1.03 0.98   5.64   9.05 1.00     2071     1664
##    y1[20]   9.94   9.94 1.04 0.99   8.27  11.66 1.00     1898     1952
##    y1[21]   9.61   9.58 1.04 1.06   7.94  11.36 1.00     2022     1927
##    y1[22]   0.34   0.30 1.11 1.10  -1.41   2.20 1.00     1528     1460
##    y1[23]   7.39   7.40 1.02 0.99   5.70   9.00 1.00     1846     1958
##    y1[24]   3.12   3.10 1.05 1.03   1.41   4.82 1.00     2090     1968
##    y1[25]   4.47   4.47 1.06 1.08   2.72   6.21 1.00     2181     1970
##    y1[26]   9.23   9.23 1.03 1.00   7.58  10.90 1.00     1965     1967
##    y1[27]  10.19  10.19 1.02 1.00   8.60  11.88 1.00     1957     1735
##    y1[28]   1.09   1.06 1.04 1.03  -0.62   2.86 1.00     1559     1734
##    y1[29]  10.39  10.40 1.06 1.04   8.64  12.12 1.00     1956     1855
##    y1[30]   5.90   5.89 1.03 1.03   4.25   7.60 1.00     1536     1648
##    m1[1]    5.15   5.13 0.49 0.50   4.37   6.00 1.00     1114     1261
##    m1[2]    4.69   4.69 0.32 0.31   4.17   5.21 1.00     1491     1259
##    m1[3]    7.16   7.15 0.44 0.43   6.43   7.89 1.00      972     1186
##    m1[4]    6.19   6.17 0.52 0.53   5.37   7.08 1.00     1195     1307
##    m1[5]    3.17   3.18 0.33 0.33   2.62   3.70 1.00      760     1026
##    m1[6]    5.41   5.41 0.36 0.35   4.81   6.01 1.00      849     1222
##    m1[7]    0.93   0.94 0.39 0.38   0.29   1.55 1.00      835      937
##    m1[8]    5.44   5.44 0.37 0.35   4.84   6.05 1.00      851     1224
##    m1[9]    6.16   6.15 0.26 0.25   5.74   6.57 1.00     1646     1334
##    m1[10]   2.64   2.66 0.34 0.33   2.09   3.19 1.00      766      977
##    m1[11]   8.84   8.85 0.27 0.26   8.39   9.29 1.00     1709     1267
##    m1[12]   8.12   8.12 0.25 0.24   7.68   8.53 1.00     1805     1324
##    m1[13]   3.83   3.81 0.48 0.48   3.09   4.65 1.00     1016     1194
##    m1[14]   9.72   9.73 0.31 0.29   9.20  10.24 1.00     1581     1274
##    m1[15]   3.96   3.96 0.36 0.36   3.38   4.56 1.00     1453     1068
##    m1[16]   5.96   5.96 0.26 0.26   5.54   6.39 1.00     1623     1373
##    m1[17]   0.59   0.60 0.40 0.39  -0.07   1.24 1.00      855      936
##    m1[18]   7.31   7.31 0.24 0.23   6.91   7.72 1.00     1800     1664
##    m1[19]   7.37   7.37 0.24 0.23   6.95   7.77 1.00     1805     1608
##    m1[20]   9.91   9.92 0.32 0.31   9.37  10.44 1.00     1558     1269
##    m1[21]   9.62   9.63 0.31 0.28   9.10  10.12 1.00     1596     1291
##    m1[22]   0.29   0.27 0.55 0.56  -0.58   1.21 1.01      905      951
##    m1[23]   7.38   7.38 0.24 0.23   6.97   7.79 1.00     1806     1608
##    m1[24]   3.12   3.12 0.41 0.41   2.44   3.80 1.00     1402     1028
##    m1[25]   4.48   4.48 0.33 0.33   3.94   5.01 1.00     1494     1195
##    m1[26]   9.21   9.22 0.29 0.26   8.73   9.68 1.00     1650     1287
##    m1[27]  10.17  10.17 0.34 0.32   9.60  10.72 1.00     1533     1297
##    m1[28]   1.11   1.13 0.38 0.37   0.49   1.72 1.00      825      928
##    m1[29]  10.39  10.39 0.35 0.33   9.80  10.96 1.00     1511     1365
##    m1[30]   5.87   5.87 0.38 0.38   5.23   6.50 1.00      880     1361

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)

plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)

qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')+
  geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')

qplot(data=tb,1:n,y,col=c)+
  geom_point(aes(x=1:n,y=y1),tb,col='black')+
  geom_line(aes(x=1:n,y=l5),tb,col='gray')+
  geom_line(aes(x=1:n,y=u5),tb,col='gray')

ex6

generalized linear regression

poisson regression

objective variable [0,Infinity), integer. variance of error is near to mean  
(normal linear regression, correlation is near to 1,-1, variance of error is constant)  

# of samples is N,  
log li=sum(bj*xji),j=[0,K],i=[1,N]  
yi~Po(li)  
 or  
li=sum(bj xji),j=[0,k]  
yi~Po(exp li)  

when xj -> xj +1, y -> y* exp bj   

ex6-1.stan

data{
  int N;
  int K;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] l=X*b;
  y~poisson_log(l);
}
generated quantities{
  array[N] int y1;
  vector[N] l1=X*b;
  for(i in 1:N){
    y1[i]=poisson_log_rng(l1[i]);
  }
}
n=30
tb=tibble(x=runif(n,-1,1),c=sample(c('a','b'),n,replace=T),
          y=rpois(n,exp(x+(c=='b')*0.5)))
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

glm(f,tb,family='poisson')
## 
## Call:  glm(formula = f, family = "poisson", data = tb)
## 
## Coefficients:
## (Intercept)            x           cb  
##       0.347        1.183        0.197  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  27 Residual
## Null Deviance:       52.2 
## Residual Deviance: 32    AIC: 96.5
X=model.matrix(f,tb)

data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -71.4455 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       10      -12.1764      0.001982   0.000665165           1           1       13    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -12.18
##    b[1]       0.35
##    b[2]       1.18
##    b[3]       0.20
##    y1[1]      0.00
##    y1[2]      1.00
##    y1[3]      0.00
##    y1[4]      0.00
##    y1[5]      5.00
##    y1[6]      2.00
##    y1[7]      2.00
##    y1[8]      1.00
##    y1[9]      2.00
##    y1[10]     0.00
##    y1[11]     3.00
##    y1[12]     6.00
##    y1[13]     2.00
##    y1[14]     0.00
##    y1[15]     7.00
##    y1[16]     1.00
##    y1[17]     4.00
##    y1[18]     0.00
##    y1[19]     7.00
##    y1[20]     6.00
##    y1[21]     1.00
##    y1[22]     1.00
##    y1[23]     2.00
##    y1[24]     0.00
##    y1[25]     3.00
##    y1[26]     2.00
##    y1[27]     1.00
##    y1[28]     1.00
##    y1[29]     2.00
##    y1[30]     0.00
##    l1[1]     -0.07
##    l1[2]     -0.58
##    l1[3]      0.30
##    l1[4]     -0.06
##    l1[5]      1.66
##    l1[6]      0.31
##    l1[7]      0.78
##    l1[8]     -0.15
##    l1[9]      0.38
##    l1[10]    -0.15
##    l1[11]     1.45
##    l1[12]     1.19
##    l1[13]    -0.24
##    l1[14]    -0.39
##    l1[15]     0.86
##    l1[16]    -0.16
##    l1[17]     0.75
##    l1[18]    -0.20
##    l1[19]     1.12
##    l1[20]     1.29
##    l1[21]     0.33
##    l1[22]     0.31
##    l1[23]     1.25
##    l1[24]    -0.58
##    l1[25]     0.86
##    l1[26]     0.53
##    l1[27]     0.41
##    l1[28]     0.29
##    l1[29]     0.48
##    l1[30]     0.07
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='l1',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -13.65 -13.36 1.20 1.00 -15.91 -12.35 1.01      718     1031
##    b[1]     0.32   0.32 0.25 0.25  -0.10   0.72 1.00      967     1092
##    b[2]     1.20   1.20 0.27 0.27   0.76   1.65 1.00     1267     1023
##    b[3]     0.20   0.20 0.29 0.29  -0.27   0.67 1.00      951     1118
##    y1[1]    0.95   1.00 1.01 1.48   0.00   3.00 1.00     2049     2005
##    y1[2]    0.55   0.00 0.75 0.00   0.00   2.00 1.00     2024     1959
##    y1[3]    1.36   1.00 1.22 1.48   0.00   4.00 1.00     1819     1840
##    y1[4]    0.95   1.00 1.00 1.48   0.00   3.00 1.00     1967     1990
##    y1[5]    5.24   5.00 2.61 2.97   1.00  10.00 1.00     1873     1910
##    y1[6]    1.38   1.00 1.24 1.48   0.00   4.00 1.00     1908     1816
##    y1[7]    2.16   2.00 1.60 1.48   0.00   5.00 1.00     1748     1786
##    y1[8]    0.90   1.00 0.99 1.48   0.00   3.00 1.00     2033     1940
##    y1[9]    1.46   1.00 1.25 1.48   0.00   4.00 1.00     1842     1825
##    y1[10]   0.87   1.00 0.97 1.48   0.00   3.00 1.00     2077     2064
##    y1[11]   4.34   4.00 2.50 2.97   1.00   9.00 1.00     1608     1335
##    y1[12]   3.25   3.00 1.85 1.48   1.00   7.00 1.00     2189     1996
##    y1[13]   0.80   1.00 0.94 1.48   0.00   3.00 1.00     1931     1836
##    y1[14]   0.68   0.00 0.86 0.00   0.00   2.00 1.00     1800     1768
##    y1[15]   2.40   2.00 1.57 1.48   0.00   5.00 1.00     2062     1914
##    y1[16]   0.88   1.00 0.97 1.48   0.00   3.00 1.00     1846     1716
##    y1[17]   2.08   2.00 1.50 1.48   0.00   5.00 1.00     2095     1971
##    y1[18]   0.83   1.00 0.92 1.48   0.00   3.00 1.00     1997     2050
##    y1[19]   3.03   3.00 1.84 1.48   0.00   6.00 1.00     2056     2046
##    y1[20]   3.70   3.00 1.97 1.48   1.00   7.00 1.00     1904     1896
##    y1[21]   1.41   1.00 1.23 1.48   0.00   4.00 1.00     1890     1819
##    y1[22]   1.33   1.00 1.18 1.48   0.00   4.00 1.00     2097     1563
##    y1[23]   3.49   3.00 1.97 1.48   1.00   7.00 1.00     1881     2004
##    y1[24]   0.56   0.00 0.78 0.00   0.00   2.00 1.00     1860     1821
##    y1[25]   2.37   2.00 1.60 1.48   0.00   5.00 1.00     1797     1707
##    y1[26]   1.65   1.00 1.34 1.48   0.00   4.00 1.00     1883     1938
##    y1[27]   1.45   1.00 1.25 1.48   0.00   4.00 1.00     1945     1912
##    y1[28]   1.32   1.00 1.19 1.48   0.00   4.00 1.00     1670     1746
##    y1[29]   1.65   1.50 1.31 0.74   0.00   4.00 1.00     1917     1910
##    y1[30]   1.04   1.00 1.04 1.48   0.00   3.00 1.00     2111     2025
##    l1[1]   -0.10  -0.10 0.29 0.29  -0.58   0.35 1.00     1007     1220
##    l1[2]   -0.62  -0.60 0.37 0.38  -1.26  -0.03 1.00     1303     1037
##    l1[3]    0.27   0.27 0.21 0.21  -0.10   0.61 1.00     1517     1252
##    l1[4]   -0.09  -0.08 0.27 0.27  -0.55   0.34 1.00     1380     1142
##    l1[5]    1.64   1.65 0.24 0.25   1.23   2.03 1.00     1656     1309
##    l1[6]    0.28   0.28 0.25 0.26  -0.15   0.68 1.00      969     1077
##    l1[7]    0.76   0.76 0.24 0.25   0.37   1.16 1.00      977     1069
##    l1[8]   -0.18  -0.18 0.30 0.30  -0.68   0.29 1.00     1018     1196
##    l1[9]    0.35   0.35 0.25 0.25  -0.06   0.74 1.00      964     1086
##    l1[10]  -0.18  -0.17 0.29 0.29  -0.68   0.27 1.00     1358     1033
##    l1[11]   1.43   1.43 0.31 0.31   0.92   1.93 1.00     1116     1075
##    l1[12]   1.17   1.18 0.18 0.18   0.86   1.44 1.00     1999     1407
##    l1[13]  -0.27  -0.26 0.30 0.31  -0.80   0.21 1.00     1341     1033
##    l1[14]  -0.43  -0.42 0.33 0.34  -0.99   0.10 1.00     1051     1189
##    l1[15]   0.84   0.85 0.16 0.16   0.57   1.10 1.00     2017     1642
##    l1[16]  -0.19  -0.19 0.30 0.30  -0.70   0.28 1.00     1020     1166
##    l1[17]   0.73   0.73 0.16 0.16   0.45   0.99 1.00     1924     1721
##    l1[18]  -0.24  -0.23 0.30 0.30  -0.75   0.23 1.00     1348     1044
##    l1[19]   1.10   1.11 0.17 0.18   0.81   1.37 1.00     2041     1401
##    l1[20]   1.28   1.29 0.19 0.20   0.95   1.57 1.00     1921     1411
##    l1[21]   0.30   0.31 0.21 0.20  -0.05   0.64 1.00     1539     1280
##    l1[22]   0.28   0.28 0.21 0.21  -0.08   0.62 1.00     1522     1252
##    l1[23]   1.23   1.24 0.18 0.19   0.92   1.52 1.00     1954     1432
##    l1[24]  -0.62  -0.61 0.36 0.36  -1.24  -0.04 1.00     1075     1158
##    l1[25]   0.83   0.83 0.25 0.25   0.44   1.24 1.00      990     1101
##    l1[26]   0.51   0.51 0.18 0.18   0.20   0.80 1.00     1693     1556
##    l1[27]   0.38   0.38 0.24 0.25  -0.03   0.78 1.00      962     1145
##    l1[28]   0.26   0.26 0.21 0.21  -0.11   0.60 1.00     1512     1252
##    l1[29]   0.46   0.46 0.19 0.19   0.14   0.76 1.00     1651     1536
##    l1[30]   0.04   0.05 0.25 0.25  -0.38   0.43 1.00     1417     1174

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##     0 1 1.5 2 3 4 5
##   0 1 7   0 0 0 0 0
##   1 1 3   0 1 1 0 0
##   2 1 5   0 2 1 0 0
##   3 0 1   1 1 0 0 0
##   4 0 0   0 0 2 0 0
##   6 0 0   0 0 0 1 0
##   7 0 0   0 0 0 0 1
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##     0 1 1.5 2 3 4 5
##   0 0 2   0 0 0 0 0
##   1 1 3   0 1 0 0 0
##   2 1 1   0 1 0 0 0
##   3 0 0   0 0 0 0 0
##   4 0 0   0 0 0 0 0
##   6 0 0   0 0 0 1 0
##   7 0 0   0 0 0 0 0
## 
## , ,  = b
## 
##    
##     0 1 1.5 2 3 4 5
##   0 1 5   0 0 0 0 0
##   1 0 0   0 0 1 0 0
##   2 0 4   0 1 1 0 0
##   3 0 1   1 1 0 0 0
##   4 0 0   0 0 2 0 0
##   6 0 0   0 0 0 0 0
##   7 0 0   0 0 0 0 1
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##     0 1 2 3 5
##   0 4 4 0 0 0
##   1 4 1 1 0 0
##   2 1 6 1 1 0
##   3 1 2 0 0 0
##   4 0 0 0 2 0
##   6 0 0 0 1 0
##   7 0 0 0 0 1
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##     0 1 2 3 5
##   0 0 2 0 0 0
##   1 4 1 0 0 0
##   2 1 2 0 0 0
##   3 0 0 0 0 0
##   4 0 0 0 0 0
##   6 0 0 0 1 0
##   7 0 0 0 0 0
## 
## , ,  = b
## 
##    map
##     0 1 2 3 5
##   0 4 2 0 0 0
##   1 0 0 1 0 0
##   2 0 4 1 1 0
##   3 1 2 0 0 0
##   4 0 0 0 2 0
##   6 0 0 0 0 0
##   7 0 0 0 0 1
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

logistic regression   

# of samples is N,  
objective variable 0/1 binary  
  
probability of incident pi[0,1]  
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)  

yi~Ber(pi), 0/1 binary  

odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident  
odds ratio(x0,x1)=odds(x1)/odd(x0)  
  
when xj -> xj +1, odds ratio -> odds ratio *exp bj  

ex6-2.stan

data{
  int N;
  int K;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] z=X*b;
  y~bernoulli_logit(z);
}
generated quantities{
  array[N] int y1;
  vector[N] z1=X*b;
  for(i in 1:N){
    y1[i]=bernoulli_rng(inv_logit(z1[i]));
  }
}
n=30
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,1,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)

f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

glm(f,tb,family='binomial') # it can caluculte when all trials are once
## 
## Call:  glm(formula = f, family = "binomial", data = tb)
## 
## Coefficients:
## (Intercept)            x           cb  
##       0.262        1.222        1.483  
## 
## Degrees of Freedom: 29 Total (i.e. Null);  27 Residual
## Null Deviance:       36.7 
## Residual Deviance: 32.6  AIC: 38.6
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mdl=cmdstan_model('./ex6-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -19.4771 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        9      -16.3142   0.000328922   0.000143539       0.954       0.954       12    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -16.31
##    b[1]       0.26
##    b[2]       1.22
##    b[3]       1.48
##    y1[1]      1.00
##    y1[2]      0.00
##    y1[3]      1.00
##    y1[4]      1.00
##    y1[5]      1.00
##    y1[6]      1.00
##    y1[7]      1.00
##    y1[8]      0.00
##    y1[9]      1.00
##    y1[10]     0.00
##    y1[11]     0.00
##    y1[12]     1.00
##    y1[13]     1.00
##    y1[14]     1.00
##    y1[15]     1.00
##    y1[16]     1.00
##    y1[17]     1.00
##    y1[18]     1.00
##    y1[19]     1.00
##    y1[20]     1.00
##    y1[21]     0.00
##    y1[22]     0.00
##    y1[23]     0.00
##    y1[24]     0.00
##    y1[25]     1.00
##    y1[26]     1.00
##    y1[27]     0.00
##    y1[28]     1.00
##    y1[29]     1.00
##    y1[30]     1.00
##    z1[1]      0.47
##    z1[2]      0.06
##    z1[3]      0.02
##    z1[4]      0.83
##    z1[5]      1.87
##    z1[6]      1.14
##    z1[7]      1.46
##    z1[8]      1.09
##    z1[9]      1.29
##    z1[10]     0.00
##    z1[11]     0.02
##    z1[12]     2.49
##    z1[13]     2.36
##    z1[14]     1.19
##    z1[15]     1.12
##    z1[16]     1.86
##    z1[17]     1.13
##    z1[18]     1.36
##    z1[19]    -0.62
##    z1[20]     2.70
##    z1[21]    -0.73
##    z1[22]     0.57
##    z1[23]     0.71
##    z1[24]     1.14
##    z1[25]     1.01
##    z1[26]     0.82
##    z1[27]    -0.47
##    z1[28]     1.43
##    z1[29]     0.95
##    z1[30]     2.18
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='z1',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -17.92 -17.62 1.24 1.05 -20.35 -16.52 1.01      790     1133
##    b[1]     0.27   0.27 0.56 0.56  -0.65   1.17 1.00     1536     1568
##    b[2]     1.40   1.34 0.92 0.86  -0.04   3.04 1.00     1006     1123
##    b[3]     1.76   1.71 1.08 1.08   0.12   3.62 1.01      982     1077
##    y1[1]    0.63   1.00 0.48 0.00   0.00   1.00 1.00     1915       NA
##    y1[2]    0.53   1.00 0.50 0.00   0.00   1.00 1.00     2065       NA
##    y1[3]    0.51   1.00 0.50 0.00   0.00   1.00 1.00     1867       NA
##    y1[4]    0.71   1.00 0.45 0.00   0.00   1.00 1.00     2120       NA
##    y1[5]    0.87   1.00 0.33 0.00   0.00   1.00 1.00     1874       NA
##    y1[6]    0.76   1.00 0.42 0.00   0.00   1.00 1.00     1783       NA
##    y1[7]    0.82   1.00 0.38 0.00   0.00   1.00 1.00     2072       NA
##    y1[8]    0.75   1.00 0.43 0.00   0.00   1.00 1.00     1966       NA
##    y1[9]    0.78   1.00 0.41 0.00   0.00   1.00 1.00     1935       NA
##    y1[10]   0.47   0.00 0.50 0.00   0.00   1.00 1.00     1918       NA
##    y1[11]   0.50   1.00 0.50 0.00   0.00   1.00 1.00     1995       NA
##    y1[12]   0.92   1.00 0.27 0.00   0.00   1.00 1.00     2025       NA
##    y1[13]   0.90   1.00 0.30 0.00   0.00   1.00 1.00     1889       NA
##    y1[14]   0.77   1.00 0.42 0.00   0.00   1.00 1.00     1900       NA
##    y1[15]   0.76   1.00 0.43 0.00   0.00   1.00 1.00     1811       NA
##    y1[16]   0.87   1.00 0.34 0.00   0.00   1.00 1.00     1896       NA
##    y1[17]   0.76   1.00 0.43 0.00   0.00   1.00 1.00     1923       NA
##    y1[18]   0.79   1.00 0.41 0.00   0.00   1.00 1.00     1930       NA
##    y1[19]   0.35   0.00 0.48 0.00   0.00   1.00 1.00     1540       NA
##    y1[20]   0.93   1.00 0.26 0.00   0.00   1.00 1.00     1914       NA
##    y1[21]   0.33   0.00 0.47 0.00   0.00   1.00 1.00     1892       NA
##    y1[22]   0.62   1.00 0.49 0.00   0.00   1.00 1.00     1911       NA
##    y1[23]   0.67   1.00 0.47 0.00   0.00   1.00 1.00     2071       NA
##    y1[24]   0.76   1.00 0.43 0.00   0.00   1.00 1.00     1944       NA
##    y1[25]   0.73   1.00 0.44 0.00   0.00   1.00 1.00     2108       NA
##    y1[26]   0.71   1.00 0.45 0.00   0.00   1.00 1.00     1882       NA
##    y1[27]   0.38   0.00 0.49 0.00   0.00   1.00 1.00     1834       NA
##    y1[28]   0.82   1.00 0.38 0.00   0.00   1.00 1.00     2027       NA
##    y1[29]   0.73   1.00 0.45 0.00   0.00   1.00 1.00     1847       NA
##    y1[30]   0.90   1.00 0.30 0.00   0.00   1.00 1.00     2012       NA
##    z1[1]    0.51   0.51 0.55 0.55  -0.38   1.41 1.00     1704     1665
##    z1[2]    0.03   0.03 0.61 0.60  -0.98   1.06 1.00     1356     1567
##    z1[3]    0.00  -0.01 0.62 0.60  -1.03   1.03 1.00     1329     1542
##    z1[4]    0.92   0.90 0.63 0.64  -0.08   1.97 1.00     1727     1401
##    z1[5]    2.17   2.08 0.92 0.93   0.80   3.81 1.01     1229     1202
##    z1[6]    1.28   1.24 0.77 0.75   0.09   2.59 1.00     1608     1219
##    z1[7]    1.70   1.63 0.88 0.88   0.39   3.26 1.00     1380     1188
##    z1[8]    1.22   1.19 0.75 0.73   0.08   2.49 1.00     1628     1219
##    z1[9]    1.51   1.44 0.89 0.87   0.18   3.08 1.00     1425     1252
##    z1[10]  -0.03  -0.04 0.63 0.62  -1.09   1.02 1.00     1311     1494
##    z1[11]  -0.01  -0.01 0.62 0.60  -1.04   1.03 1.00     1328     1518
##    z1[12]   2.88   2.79 1.15 1.16   1.15   4.88 1.01      991     1170
##    z1[13]   2.73   2.65 1.09 1.11   1.09   4.67 1.01     1014     1216
##    z1[14]   1.40   1.34 0.91 0.90   0.04   2.98 1.00     1440     1251
##    z1[15]   1.32   1.27 0.92 0.93  -0.06   2.93 1.00     1451     1303
##    z1[16]   2.16   2.07 0.92 0.93   0.78   3.78 1.01     1229     1202
##    z1[17]   1.26   1.22 0.76 0.74   0.08   2.56 1.00     1616     1219
##    z1[18]   1.53   1.48 0.90 0.87   0.17   3.10 1.00     1478     1222
##    z1[19]  -0.74  -0.73 0.95 0.89  -2.35   0.81 1.00     1027     1238
##    z1[20]   3.12   3.01 1.26 1.28   1.22   5.30 1.01      966     1089
##    z1[21]  -0.86  -0.85 1.02 0.95  -2.63   0.81 1.00     1002     1238
##    z1[22]   0.63   0.62 0.56 0.56  -0.29   1.53 1.00     1743     1542
##    z1[23]   0.85   0.80 1.06 1.07  -0.81   2.68 1.00     1396     1380
##    z1[24]   1.27   1.23 0.77 0.74   0.09   2.58 1.00     1609     1219
##    z1[25]   1.12   1.10 0.70 0.70   0.01   2.30 1.00     1662     1305
##    z1[26]   0.91   0.89 0.63 0.64  -0.08   1.95 1.00     1730     1401
##    z1[27]  -0.57  -0.57 0.86 0.82  -2.04   0.86 1.00     1073     1307
##    z1[28]   1.68   1.61 0.88 0.87   0.36   3.22 1.00     1387     1188
##    z1[29]   1.06   1.03 0.68 0.68  -0.01   2.19 1.00     1684     1286
##    z1[30]   2.53   2.44 1.02 1.04   1.00   4.31 1.01     1060     1246

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##      0  1
##   0  3  6
##   1  1 20
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##      0  1
##   0  3  4
##   1  1 10
## 
## , ,  = b
## 
##    
##      0  1
##   0  0  2
##   1  0 10
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##      0  1
##   0  3  6
##   1  1 20
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##      0  1
##   0  3  4
##   1  1 10
## 
## , ,  = b
## 
##    map
##      0  1
##   0  0  2
##   1  0 10
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

multi logistic regression  

# of samples is N,  
# of trials of a sample i is mi,  
objective variable [0,n], integer  
  
probability of incident pi[0,1]  
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)  

yi~B(mi, pi), # of acutual incident  

odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident  
odds ratio(x0,x1)=odds(x1)/odd(x0)  
  
when xj -> xj +1, odds ratio -> odds ratio *exp bj  

ex6-3.stan

data{
  int N;
  int K;
  array[N] int m;
  array[N] int y;
  matrix[N,K] X;
}
parameters{
  vector[K] b;
}
model{
  vector[N] z=X*b;
  y~binomial_logit(m,z);
}
generated quantities{
  array[N] int y1;
  vector[N] z1=X*b;
  for(i in 1:N){
    y1[i]=binomial_rng(m[i],inv_logit(z1[i]));
  }
}
n=30
m=floor(runif(n,1,10)) # trials are varying (1,10)
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,m,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)

f=formula(y~x+c)
qplot(data=tb,x,y,col=c)

X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X,m=m)

mdl=cmdstan_model('./ex6-3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -96.0634 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       10      -90.0388   0.000161926   5.07823e-05      0.9345      0.9345       12    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -90.04
##    b[1]       0.23
##    b[2]       0.60
##    b[3]       1.43
##    y1[1]      5.00
##    y1[2]      8.00
##    y1[3]      3.00
##    y1[4]      3.00
##    y1[5]      7.00
##    y1[6]      2.00
##    y1[7]      5.00
##    y1[8]      2.00
##    y1[9]      1.00
##    y1[10]     1.00
##    y1[11]     2.00
##    y1[12]     2.00
##    y1[13]     5.00
##    y1[14]     1.00
##    y1[15]     6.00
##    y1[16]     3.00
##    y1[17]     3.00
##    y1[18]     1.00
##    y1[19]     2.00
##    y1[20]     8.00
##    y1[21]     0.00
##    y1[22]     2.00
##    y1[23]     4.00
##    y1[24]     8.00
##    y1[25]     3.00
##    y1[26]     2.00
##    y1[27]     4.00
##    y1[28]     5.00
##    y1[29]     4.00
##    y1[30]     2.00
##    z1[1]      0.61
##    z1[2]      2.11
##    z1[3]      1.73
##    z1[4]      0.68
##    z1[5]      2.02
##    z1[6]      0.76
##    z1[7]      1.91
##    z1[8]      0.02
##    z1[9]      1.09
##    z1[10]     0.12
##    z1[11]    -0.18
##    z1[12]     1.47
##    z1[13]     1.57
##    z1[14]     0.21
##    z1[15]    -0.06
##    z1[16]     0.18
##    z1[17]     1.45
##    z1[18]     0.55
##    z1[19]     1.60
##    z1[20]     2.10
##    z1[21]     0.50
##    z1[22]    -0.21
##    z1[23]     1.61
##    z1[24]     2.26
##    z1[25]     1.45
##    z1[26]     0.51
##    z1[27]     0.76
##    z1[28]     2.01
##    z1[29]     0.39
##    z1[30]     1.65
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,exc='z1',ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -91.57 -91.27 1.21 1.00 -93.91 -90.23 1.00      901     1189
##    b[1]     0.24   0.24 0.23 0.24  -0.13   0.62 1.00     1438     1403
##    b[2]     0.59   0.60 0.34 0.35   0.04   1.12 1.00     1491     1325
##    b[3]     1.49   1.47 0.41 0.41   0.83   2.16 1.00     1164      910
##    y1[1]    5.83   6.00 1.52 1.48   3.00   8.00 1.00     1875     1886
##    y1[2]    7.12   7.00 0.93 1.48   5.00   8.00 1.00     2016       NA
##    y1[3]    3.42   4.00 0.71 0.00   2.00   4.00 1.00     2021       NA
##    y1[4]    2.63   3.00 0.98 1.48   1.00   4.00 1.00     1885       NA
##    y1[5]    7.07   7.00 0.95 1.48   5.00   8.00 1.00     1920       NA
##    y1[6]    2.72   3.00 0.97 1.48   1.00   4.00 1.00     1969       NA
##    y1[7]    5.22   5.00 0.84 1.48   4.00   6.00 1.00     1990       NA
##    y1[8]    3.50   4.00 1.41 1.48   1.00   6.00 1.00     1867     1909
##    y1[9]    1.50   2.00 0.63 0.00   0.00   2.00 1.00     1898       NA
##    y1[10]   1.06   1.00 0.71 1.48   0.00   2.00 1.00     1934       NA
##    y1[11]   1.36   1.00 0.89 1.48   0.00   3.00 1.00     1914       NA
##    y1[12]   3.29   3.00 0.79 1.48   2.00   4.00 1.00     1811       NA
##    y1[13]   6.66   7.00 1.13 1.48   5.00   8.00 1.00     1811       NA
##    y1[14]   1.65   2.00 0.88 1.48   0.00   3.00 1.00     1992       NA
##    y1[15]   4.37   4.00 1.65 1.48   2.00   7.00 1.00     1851     1938
##    y1[16]   3.29   3.00 1.22 1.48   1.00   5.00 1.00     2077     1883
##    y1[17]   3.24   3.00 0.81 1.48   2.00   4.00 1.00     1993       NA
##    y1[18]   5.04   5.00 1.44 1.48   3.00   7.00 1.00     1994     1945
##    y1[19]   1.66   2.00 0.53 0.00   1.00   2.00 1.00     1850       NA
##    y1[20]   7.12   7.00 0.92 1.48   5.00   8.00 1.00     1931       NA
##    y1[21]   2.50   3.00 1.00 1.48   1.00   4.00 1.00     1629       NA
##    y1[22]   3.58   3.50 1.58 2.22   1.00   6.00 1.00     1647     1670
##    y1[23]   4.21   4.00 0.85 1.48   3.00   5.00 1.00     1838       NA
##    y1[24]   7.24   7.00 0.86 1.48   6.00   8.00 1.00     2019       NA
##    y1[25]   2.47   3.00 0.68 0.00   1.00   3.00 1.00     1698       NA
##    y1[26]   3.74   4.00 1.20 1.48   2.00   6.00 1.00     1818       NA
##    y1[27]   5.42   6.00 1.44 1.48   3.00   8.00 1.00     1901       NA
##    y1[28]   4.42   5.00 0.74 0.00   3.00   5.00 1.00     1942       NA
##    y1[29]   4.13   4.00 1.35 1.48   2.00   6.00 1.00     1968     1959
##    y1[30]   2.53   3.00 0.65 0.00   1.00   3.00 1.00     2003       NA
##    z1[1]    0.61   0.61 0.28 0.29   0.15   1.11 1.00     1654     1585
##    z1[2]    2.17   2.15 0.38 0.38   1.57   2.82 1.00     1633     1235
##    z1[3]    1.79   1.77 0.34 0.33   1.27   2.37 1.00     1453     1346
##    z1[4]    0.68   0.68 0.31 0.31   0.17   1.20 1.00     1643     1562
##    z1[5]    2.07   2.06 0.36 0.35   1.52   2.68 1.00     1601     1285
##    z1[6]    0.76   0.76 0.34 0.35   0.20   1.34 1.00     1635     1563
##    z1[7]    1.97   1.95 0.35 0.34   1.43   2.55 1.00     1555     1362
##    z1[8]    0.03   0.03 0.28 0.28  -0.43   0.49 1.00     1499     1506
##    z1[9]    1.15   1.15 0.52 0.53   0.31   2.02 1.00     1369     1388
##    z1[10]   0.13   0.13 0.25 0.26  -0.28   0.54 1.00     1461     1416
##    z1[11]  -0.17  -0.17 0.35 0.36  -0.76   0.42 1.00     1549     1275
##    z1[12]   1.53   1.52 0.38 0.38   0.93   2.16 1.00     1367     1400
##    z1[13]   1.63   1.62 0.36 0.36   1.07   2.24 1.00     1392     1346
##    z1[14]   0.21   0.21 0.23 0.24  -0.16   0.60 1.00     1442     1393
##    z1[15]  -0.05  -0.06 0.31 0.31  -0.57   0.46 1.00     1547     1407
##    z1[16]   0.18   0.18 0.24 0.25  -0.20   0.57 1.00     1442     1568
##    z1[17]   1.51   1.50 0.39 0.39   0.89   2.16 1.00     1366     1444
##    z1[18]   0.55   0.55 0.26 0.26   0.12   1.01 1.00     1639     1535
##    z1[19]   1.65   1.64 0.36 0.36   1.10   2.25 1.00     1398     1346
##    z1[20]   2.15   2.14 0.38 0.38   1.56   2.80 1.00     1631     1251
##    z1[21]   0.50   0.49 0.25 0.25   0.09   0.94 1.00     1619     1514
##    z1[22]  -0.20  -0.20 0.37 0.38  -0.81   0.41 1.00     1546     1275
##    z1[23]   1.67   1.65 0.35 0.35   1.12   2.27 1.00     1403     1345
##    z1[24]   2.31   2.30 0.43 0.43   1.63   3.05 1.00     1643     1303
##    z1[25]   1.51   1.49 0.39 0.39   0.89   2.16 1.00     1365     1422
##    z1[26]   0.51   0.51 0.25 0.25   0.10   0.96 1.00     1625     1562
##    z1[27]   0.76   0.76 0.34 0.35   0.20   1.34 1.00     1636     1589
##    z1[28]   2.06   2.05 0.36 0.35   1.51   2.66 1.00     1598     1235
##    z1[29]   0.39   0.39 0.23 0.23   0.02   0.79 1.00     1531     1369
##    z1[30]   1.71   1.69 0.35 0.34   1.17   2.30 1.00     1419     1326

y1=mcmc$draws('y1')
smy=summary(y1)

table(tb$y,smy$median)
##    
##     1 2 3 3.5 4 5 6 7
##   0 0 1 0   0 0 0 0 0
##   1 2 0 0   0 0 0 0 0
##   2 0 1 2   0 2 0 0 0
##   3 0 1 4   0 2 0 0 0
##   4 0 0 2   0 1 1 0 0
##   5 0 0 0   0 1 2 1 1
##   6 0 0 0   1 0 0 0 0
##   7 0 0 0   0 0 0 1 1
##   8 0 0 0   0 0 0 0 3
cat('\n')
table(tb$y,smy$median,tb$c)
## , ,  = a
## 
##    
##     1 2 3 3.5 4 5 6 7
##   0 0 0 0   0 0 0 0 0
##   1 2 0 0   0 0 0 0 0
##   2 0 0 2   0 1 0 0 0
##   3 0 1 1   0 2 0 0 0
##   4 0 0 1   0 0 1 0 0
##   5 0 0 0   0 1 0 1 0
##   6 0 0 0   1 0 0 0 0
##   7 0 0 0   0 0 0 1 0
##   8 0 0 0   0 0 0 0 0
## 
## , ,  = b
## 
##    
##     1 2 3 3.5 4 5 6 7
##   0 0 1 0   0 0 0 0 0
##   1 0 0 0   0 0 0 0 0
##   2 0 1 0   0 1 0 0 0
##   3 0 0 3   0 0 0 0 0
##   4 0 0 1   0 1 0 0 0
##   5 0 0 0   0 0 2 0 1
##   6 0 0 0   0 0 0 0 0
##   7 0 0 0   0 0 0 0 1
##   8 0 0 0   0 0 0 0 3
par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

map=c()
for(i in 1:n){
  a=table(y1[,,i])
  map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)
##    map
##     1 2 3 4 5 6 7 8
##   0 0 1 0 0 0 0 0 0
##   1 2 0 0 0 0 0 0 0
##   2 0 1 2 1 1 0 0 0
##   3 0 1 3 3 0 0 0 0
##   4 0 0 1 2 1 0 0 0
##   5 0 0 0 1 1 2 1 0
##   6 0 0 1 0 0 0 0 0
##   7 0 0 0 0 0 1 0 1
##   8 0 0 0 0 0 0 0 3
cat('\n')
table(tb$y,map,tb$c)
## , ,  = a
## 
##    map
##     1 2 3 4 5 6 7 8
##   0 0 0 0 0 0 0 0 0
##   1 2 0 0 0 0 0 0 0
##   2 0 0 2 1 0 0 0 0
##   3 0 1 1 2 0 0 0 0
##   4 0 0 1 0 1 0 0 0
##   5 0 0 0 1 0 1 0 0
##   6 0 0 1 0 0 0 0 0
##   7 0 0 0 0 0 1 0 0
##   8 0 0 0 0 0 0 0 0
## 
## , ,  = b
## 
##    map
##     1 2 3 4 5 6 7 8
##   0 0 1 0 0 0 0 0 0
##   1 0 0 0 0 0 0 0 0
##   2 0 1 0 0 1 0 0 0
##   3 0 0 2 1 0 0 0 0
##   4 0 0 0 2 0 0 0 0
##   5 0 0 0 0 1 1 1 0
##   6 0 0 0 0 0 0 0 0
##   7 0 0 0 0 0 0 0 1
##   8 0 0 0 0 0 0 0 3
par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')

ex7 interaction of variables

n=50
mdl=cmdstan_model('./ex5.stan')

tb=tibble(x=runif(n,-3,3),
          ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
          y=rnorm(n,x+(ca=='b')+(cb=='b')-(ca=='b')*(cb=='b'),1))

grid.arrange(qplot(data=tb,x,y,col=ca),
             qplot(data=tb,x,y,col=cb),ncol=2)

fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
mle

mcmc=goMCMC(mdl,data)

seeMCMC(mcmc,exc='m',ch=F)

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

grid.arrange(
  qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
grid.arrange(
  qplot(data=tb,1:n,y,col=ca)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  qplot(data=tb,1:n,y,col=cb)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  ncol=2
)
}


f0=formula(y~x+ca)
f1=formula(y~x+ca+cb)
f2=formula(y~x+ca*cb)

fn(f0)
## Initial log joint probability = -60.0659 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       15      -20.7022   0.000118055    0.00119508           1           1       19    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -22.83 -22.56 1.44 1.32 -25.52 -21.12 1.00      899     1164
##    b[1]     0.14   0.14 0.21 0.20  -0.20   0.49 1.00      607      723
##    b[2]     1.22   1.22 0.08 0.08   1.09   1.36 1.00     2429     1901
##    b[3]     1.11   1.12 0.29 0.27   0.64   1.58 1.00      493      661
##    s        0.97   0.96 0.10 0.10   0.82   1.14 1.00     1622     1486
##    y1[1]   -0.27  -0.27 0.97 0.95  -1.82   1.38 1.00     1917     1996
##    y1[2]   -2.41  -2.40 1.04 1.04  -4.13  -0.73 1.00     2205     2024
##    y1[3]   -1.25  -1.26 0.99 0.99  -2.86   0.39 1.00     2175     2100
##    y1[4]   -0.71  -0.72 1.00 0.98  -2.39   0.97 1.00     1770     1836
##    y1[5]   -1.15  -1.16 1.00 0.97  -2.77   0.52 1.00     1776     1757
##    y1[6]    2.91   2.90 0.98 0.99   1.31   4.45 1.00     1776     1555
##    y1[7]   -1.70  -1.67 0.99 1.00  -3.35  -0.11 1.00     1993     1942
##    y1[8]    2.90   2.90 0.99 0.94   1.26   4.52 1.00     2028     1933
##    y1[9]    1.77   1.76 1.00 0.97   0.14   3.43 1.00     1970     2008
##    y1[10]   0.41   0.38 0.98 0.98  -1.19   2.04 1.00     1908     1783
##    y1[11]  -0.43  -0.43 0.99 0.98  -2.02   1.26 1.00     1947     1882
##    y1[12]   2.23   2.24 0.98 0.95   0.62   3.88 1.00     2039     2013
##    y1[13]   4.49   4.51 1.06 1.03   2.71   6.18 1.00     1721     1805
##    y1[14]   1.93   1.94 1.01 0.97   0.26   3.54 1.00     1924     2010
##    y1[15]   0.72   0.70 0.99 0.95  -0.90   2.37 1.00     1775     1801
##    y1[16]   0.65   0.64 1.03 1.00  -1.03   2.37 1.00     2033     1831
##    y1[17]   3.46   3.48 1.02 0.98   1.78   5.11 1.00     1952     1724
##    y1[18]   2.63   2.61 1.00 0.99   1.00   4.30 1.00     1716     1874
##    y1[19]  -1.02  -1.04 1.00 0.97  -2.66   0.64 1.00     2048     1817
##    y1[20]  -2.20  -2.15 1.05 1.04  -3.99  -0.50 1.00     2032     1891
##    y1[21]  -3.29  -3.25 1.02 0.98  -4.98  -1.62 1.00     1811     1852
##    y1[22]   2.12   2.14 0.99 0.99   0.48   3.72 1.00     2008     1891
##    y1[23]   3.59   3.58 1.01 1.03   1.92   5.25 1.00     1853     1802
##    y1[24]   3.51   3.50 0.99 1.00   1.93   5.20 1.00     1944     1952
##    y1[25]   2.60   2.59 1.00 0.98   0.96   4.25 1.00     2017     1807
##    y1[26]   4.77   4.76 1.02 1.00   3.13   6.41 1.00     1923     1852
##    y1[27]  -1.44  -1.43 1.01 1.00  -3.08   0.15 1.00     2014     1730
##    y1[28]   1.31   1.30 1.01 1.01  -0.35   2.96 1.00     1906     1890
##    y1[29]  -1.65  -1.68 1.00 1.00  -3.33   0.02 1.00     1891     1892
##    y1[30]  -1.99  -2.00 0.98 0.97  -3.57  -0.34 1.00     1748     1790
##    y1[31]   1.33   1.35 1.03 1.03  -0.34   2.99 1.00     1690     2014
##    y1[32]  -0.21  -0.21 1.00 0.98  -1.84   1.42 1.00     1754     1593
##    y1[33]  -0.54  -0.56 0.99 1.01  -2.14   1.02 1.00     2027     1801
##    y1[34]  -0.54  -0.55 1.01 0.99  -2.23   1.14 1.00     1806     1970
##    y1[35]   1.52   1.50 1.02 1.00  -0.14   3.21 1.00     1851     1883
##    y1[36]   0.66   0.66 0.98 0.95  -0.91   2.31 1.00     1760     1909
##    y1[37]  -1.82  -1.82 0.99 0.93  -3.43  -0.18 1.00     2189     1917
##    y1[38]   0.87   0.88 1.00 1.01  -0.74   2.51 1.00     1993     1788
##    y1[39]   4.27   4.31 1.02 1.01   2.57   5.92 1.00     1866     1953
##    y1[40]   3.49   3.51 1.02 1.01   1.71   5.16 1.00     1909     1843
##    y1[41]  -0.24  -0.22 1.01 0.99  -1.95   1.41 1.00     1914     1970
##    y1[42]   2.24   2.24 1.01 1.02   0.64   3.89 1.00     1830     1857
##    y1[43]   3.43   3.43 1.01 1.03   1.81   5.08 1.00     1884     1959
##    y1[44]   1.16   1.18 1.01 0.94  -0.52   2.82 1.00     1813     1843
##    y1[45]  -0.75  -0.73 1.01 0.98  -2.43   0.95 1.00     1887     1949
##    y1[46]   2.15   2.15 0.98 0.94   0.60   3.74 1.00     1905     1843
##    y1[47]   3.60   3.60 1.03 1.06   1.94   5.31 1.00     1949     1939
##    y1[48]  -1.70  -1.69 0.99 0.96  -3.35  -0.07 1.00     1917     1821
##    y1[49]  -0.17  -0.18 0.99 0.97  -1.83   1.41 1.00     2093     1969
##    y1[50]   2.31   2.30 1.02 1.05   0.64   3.98 1.00     1756     1681
##    m1[1]   -0.29  -0.29 0.22 0.21  -0.65   0.07 1.00      595      728
##    m1[2]   -2.40  -2.40 0.29 0.28  -2.90  -1.92 1.00     2403     1543
##    m1[3]   -1.25  -1.25 0.24 0.23  -1.65  -0.86 1.00     1973     1616
##    m1[4]   -0.72  -0.72 0.23 0.22  -1.10  -0.34 1.00      600      802
##    m1[5]   -1.16  -1.16 0.25 0.23  -1.57  -0.76 1.00      618      829
##    m1[6]    2.92   2.92 0.23 0.22   2.55   3.30 1.00     1023     1382
##    m1[7]   -1.69  -1.68 0.26 0.25  -2.12  -1.27 1.00     2161     1618
##    m1[8]    2.89   2.88 0.22 0.22   2.52   3.26 1.00     1020     1382
##    m1[9]    1.79   1.78 0.22 0.22   1.42   2.15 1.00      858     1027
##    m1[10]   0.41   0.41 0.21 0.20   0.08   0.75 1.00      622      703
##    m1[11]  -0.45  -0.45 0.23 0.22  -0.81  -0.08 1.00      595      776
##    m1[12]   2.22   2.21 0.23 0.23   1.84   2.60 1.00      999     1099
##    m1[13]   4.45   4.45 0.30 0.29   3.96   4.94 1.00     1208     1554
##    m1[14]   1.94   1.93 0.20 0.19   1.61   2.26 1.00      983     1171
##    m1[15]   0.68   0.67 0.21 0.20   0.35   1.02 1.00      646      751
##    m1[16]   0.65   0.64 0.19 0.18   0.34   0.96 1.00     1185     1340
##    m1[17]   3.47   3.47 0.25 0.24   3.06   3.87 1.00     1082     1382
##    m1[18]   2.65   2.64 0.25 0.25   2.25   3.04 1.00     1173     1211
##    m1[19]  -1.03  -1.03 0.23 0.23  -1.42  -0.66 1.00     1879     1638
##    m1[20]  -2.19  -2.18 0.28 0.27  -2.67  -1.73 1.00     2341     1543
##    m1[21]  -3.29  -3.28 0.34 0.34  -3.87  -2.72 1.00      793     1095
##    m1[22]   2.12   2.11 0.20 0.19   1.79   2.44 1.00      982     1283
##    m1[23]   3.60   3.60 0.25 0.25   3.18   4.02 1.00     1096     1445
##    m1[24]   3.52   3.51 0.25 0.25   3.10   3.92 1.00     1087     1402
##    m1[25]   2.59   2.59 0.24 0.24   2.19   2.98 1.00     1149     1211
##    m1[26]   4.80   4.80 0.31 0.31   4.28   5.32 1.00     1240     1597
##    m1[27]  -1.45  -1.44 0.25 0.24  -1.86  -1.05 1.00     2055     1618
##    m1[28]   1.33   1.32 0.21 0.21   0.99   1.67 1.00      746      879
##    m1[29]  -1.61  -1.60 0.25 0.25  -2.04  -1.19 1.00     2130     1618
##    m1[30]  -1.98  -1.98 0.28 0.27  -2.45  -1.53 1.00      673      979
##    m1[31]   1.32   1.31 0.19 0.18   1.00   1.63 1.00     1034     1179
##    m1[32]  -0.21  -0.21 0.22 0.21  -0.56   0.15 1.00      596      718
##    m1[33]  -0.54  -0.54 0.21 0.21  -0.90  -0.19 1.00     1667     1562
##    m1[34]  -0.49  -0.49 0.23 0.21  -0.86  -0.12 1.00      596      776
##    m1[35]   1.51   1.51 0.22 0.21   1.17   1.86 1.00      787      950
##    m1[36]   0.68   0.68 0.21 0.20   0.35   1.02 1.00      647      751
##    m1[37]  -1.85  -1.84 0.27 0.26  -2.30  -1.42 1.00     2224     1645
##    m1[38]   0.87   0.87 0.19 0.18   0.57   1.18 1.00     1122     1194
##    m1[39]   4.26   4.26 0.29 0.28   3.79   4.73 1.00     1181     1552
##    m1[40]   3.50   3.50 0.28 0.27   3.04   3.95 1.00     1555     1645
##    m1[41]  -0.22  -0.22 0.22 0.21  -0.57   0.14 1.00      595      718
##    m1[42]   2.23   2.22 0.23 0.23   1.85   2.61 1.00     1003     1099
##    m1[43]   3.42   3.42 0.28 0.27   2.97   3.87 1.00     1518     1587
##    m1[44]   1.18   1.18 0.19 0.18   0.87   1.49 1.00     1055     1210
##    m1[45]  -0.76  -0.76 0.22 0.21  -1.14  -0.40 1.00     1760     1643
##    m1[46]   2.13   2.12 0.20 0.19   1.80   2.46 1.00      982     1286
##    m1[47]   3.59   3.59 0.25 0.25   3.18   4.01 1.00     1095     1445
##    m1[48]  -1.71  -1.71 0.26 0.25  -2.15  -1.29 1.00     2169     1618
##    m1[49]  -0.19  -0.19 0.20 0.20  -0.54   0.14 1.00     1523     1465
##    m1[50]   2.32   2.31 0.24 0.23   1.94   2.70 1.00     1036     1098

fn(f1)
## Initial log joint probability = -1750.03 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       24      -19.2198   4.63272e-05   0.000150403           1           1       39    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -21.99 -21.62 1.72 1.55 -25.26 -19.91 1.01      797      989
##    b[1]    -0.07  -0.07 0.25 0.24  -0.50   0.34 1.00      844      972
##    b[2]     1.24   1.24 0.08 0.08   1.10   1.37 1.00     2191     1424
##    b[3]     1.09   1.09 0.29 0.28   0.64   1.56 1.00      858      871
##    b[4]     0.44   0.45 0.27 0.26  -0.02   0.89 1.00     1009     1013
##    s        0.95   0.94 0.10 0.10   0.81   1.14 1.00     2295     1549
##    y1[1]   -0.50  -0.51 0.98 0.96  -2.07   1.08 1.00     1965     1866
##    y1[2]   -2.65  -2.66 1.03 1.03  -4.32  -0.94 1.00     2034     1685
##    y1[3]   -1.49  -1.47 1.01 1.00  -3.14   0.14 1.00     1769     1856
##    y1[4]   -0.51  -0.51 0.97 0.97  -2.14   1.05 1.00     1775     1786
##    y1[5]   -0.96  -0.97 0.99 0.95  -2.64   0.65 1.00     1707     1761
##    y1[6]    2.67   2.69 1.01 1.01   0.99   4.34 1.00     1848     1827
##    y1[7]   -1.54  -1.55 1.01 1.02  -3.19   0.13 1.00     1932     1936
##    y1[8]    3.12   3.11 0.98 0.97   1.50   4.76 1.00     1905     1971
##    y1[9]    2.04   2.06 1.02 1.01   0.38   3.71 1.00     1879     1925
##    y1[10]   0.62   0.63 0.97 0.98  -0.94   2.20 1.00     1930     1989
##    y1[11]  -0.67  -0.69 1.00 0.99  -2.33   1.01 1.00     1878     1730
##    y1[12]   2.03   2.04 1.02 1.00   0.31   3.66 1.00     1996     2093
##    y1[13]   4.66   4.66 1.03 1.01   2.97   6.35 1.00     1959     1973
##    y1[14]   1.66   1.65 0.98 0.98   0.05   3.24 1.00     1940     1972
##    y1[15]   0.92   0.94 0.99 1.00  -0.72   2.48 1.00     1900     1890
##    y1[16]   0.85   0.86 0.97 0.99  -0.77   2.41 1.00     2094     1743
##    y1[17]   3.69   3.70 1.02 0.99   2.06   5.39 1.00     1907     1828
##    y1[18]   2.91   2.88 1.02 1.01   1.27   4.59 1.00     2080     2059
##    y1[19]  -1.27  -1.26 1.01 1.02  -2.94   0.35 1.00     1952     1945
##    y1[20]  -2.00  -1.97 1.00 0.97  -3.66  -0.41 1.00     1898     1813
##    y1[21]  -3.51  -3.53 1.02 0.99  -5.20  -1.87 1.00     1746     1783
##    y1[22]   1.84   1.84 0.99 0.99   0.19   3.47 1.00     1919     1858
##    y1[23]   3.38   3.37 1.00 0.96   1.76   5.10 1.00     1966     1852
##    y1[24]   3.29   3.31 1.01 0.98   1.63   4.93 1.00     1968     1801
##    y1[25]   2.41   2.42 1.00 0.97   0.72   4.09 1.00     1737     1931
##    y1[26]   5.09   5.12 1.03 1.00   3.39   6.76 1.00     1989     1984
##    y1[27]  -1.26  -1.28 1.01 1.00  -2.86   0.43 1.00     1955     1856
##    y1[28]   1.11   1.11 0.99 0.97  -0.57   2.72 1.00     1856     1877
##    y1[29]  -1.49  -1.48 1.00 1.01  -3.14   0.15 1.00     1907     1832
##    y1[30]  -1.75  -1.79 1.00 0.99  -3.34  -0.06 1.00     1806     1483
##    y1[31]   1.09   1.10 0.96 0.94  -0.47   2.66 1.00     1883     1907
##    y1[32]  -0.41  -0.42 0.98 0.95  -1.99   1.22 1.00     1683     1816
##    y1[33]  -0.80  -0.77 1.00 0.97  -2.48   0.84 1.00     1987     1740
##    y1[34]  -0.71  -0.73 0.99 0.98  -2.29   0.95 1.00     1876     1927
##    y1[35]   1.29   1.28 0.96 0.94  -0.23   2.90 1.00     1578     1909
##    y1[36]   0.92   0.91 1.00 0.97  -0.72   2.53 1.00     1910     1970
##    y1[37]  -1.68  -1.68 0.98 0.99  -3.35  -0.03 1.00     2019     2010
##    y1[38]   0.61   0.60 1.00 0.97  -1.03   2.22 1.00     1917     1865
##    y1[39]   4.05   4.06 1.00 0.98   2.41   5.70 1.00     2006     2084
##    y1[40]   3.83   3.84 1.02 0.99   2.12   5.48 1.00     1732     1854
##    y1[41]  -0.42  -0.45 0.99 0.94  -2.00   1.19 1.00     2013     1806
##    y1[42]   2.49   2.50 0.99 0.96   0.83   4.11 1.00     1686     1934
##    y1[43]   3.26   3.27 0.99 0.99   1.65   4.87 1.00     1807     1850
##    y1[44]   1.39   1.39 0.98 0.97  -0.20   3.04 1.00     2083     1965
##    y1[45]  -0.56  -0.56 0.97 0.96  -2.22   1.01 1.00     1930     1972
##    y1[46]   2.35   2.34 1.01 1.00   0.61   4.01 1.00     1899     1902
##    y1[47]   3.38   3.37 0.98 0.94   1.73   4.94 1.00     1682     1810
##    y1[48]  -1.53  -1.53 0.98 0.96  -3.14   0.08 1.00     1915     1815
##    y1[49]  -0.01  -0.04 0.99 0.98  -1.62   1.64 1.00     2119     1844
##    y1[50]   2.12   2.13 0.98 1.00   0.56   3.70 1.00     1900     1774
##    m1[1]   -0.51  -0.51 0.26 0.25  -0.95  -0.10 1.00      845      953
##    m1[2]   -2.68  -2.68 0.34 0.34  -3.25  -2.12 1.00     1846     1847
##    m1[3]   -1.51  -1.51 0.29 0.28  -1.99  -1.03 1.00     1613     1473
##    m1[4]   -0.50  -0.49 0.26 0.25  -0.93  -0.07 1.00      957     1180
##    m1[5]   -0.94  -0.95 0.27 0.26  -1.38  -0.50 1.00      966     1187
##    m1[6]    2.70   2.70 0.26 0.25   2.28   3.12 1.00     1119     1313
##    m1[7]   -1.51  -1.51 0.28 0.27  -1.98  -1.07 1.00     2459     1807
##    m1[8]    3.11   3.11 0.27 0.26   2.66   3.54 1.00     1568     1294
##    m1[9]    2.03   2.03 0.27 0.26   1.60   2.46 1.00     1246     1590
##    m1[10]   0.65   0.64 0.25 0.24   0.22   1.06 1.00     1024     1254
##    m1[11]  -0.67  -0.66 0.26 0.26  -1.12  -0.25 1.00      850      950
##    m1[12]   2.02   2.02 0.26 0.26   1.59   2.47 1.00     1028     1271
##    m1[13]   4.69   4.70 0.34 0.33   4.13   5.24 1.00     1637     1461
##    m1[14]   1.70   1.71 0.24 0.24   1.31   2.08 1.00     1107     1177
##    m1[15]   0.91   0.91 0.25 0.24   0.49   1.33 1.00     1054     1346
##    m1[16]   0.85   0.85 0.22 0.21   0.48   1.20 1.00     1814     1416
##    m1[17]   3.70   3.70 0.29 0.28   3.20   4.17 1.00     1584     1327
##    m1[18]   2.90   2.90 0.29 0.28   2.43   3.39 1.00     1433     1616
##    m1[19]  -1.30  -1.29 0.28 0.28  -1.75  -0.83 1.00     1569     1526
##    m1[20]  -2.02  -2.01 0.30 0.30  -2.52  -1.55 1.00     2568     1874
##    m1[21]  -3.54  -3.53 0.37 0.37  -4.16  -2.93 1.00     1077     1382
##    m1[22]   1.89   1.89 0.24 0.24   1.49   2.27 1.00     1103     1176
##    m1[23]   3.38   3.39 0.29 0.28   2.92   3.85 1.00     1151     1482
##    m1[24]   3.30   3.31 0.28 0.27   2.84   3.76 1.00     1146     1482
##    m1[25]   2.40   2.40 0.27 0.27   1.94   2.85 1.00     1094     1277
##    m1[26]   5.04   5.05 0.35 0.35   4.45   5.62 1.00     1659     1525
##    m1[27]  -1.27  -1.26 0.27 0.26  -1.71  -0.84 1.00     2406     1794
##    m1[28]   1.13   1.13 0.25 0.24   0.70   1.54 1.00      910     1100
##    m1[29]  -1.43  -1.43 0.27 0.27  -1.89  -1.00 1.00     2441     1841
##    m1[30]  -1.77  -1.77 0.30 0.29  -2.26  -1.30 1.00      999     1144
##    m1[31]   1.08   1.08 0.24 0.23   0.69   1.46 1.00     1150     1296
##    m1[32]  -0.43  -0.42 0.26 0.25  -0.86  -0.02 1.00      844      961
##    m1[33]  -0.80  -0.79 0.26 0.26  -1.23  -0.37 1.00     1474     1491
##    m1[34]  -0.71  -0.71 0.26 0.26  -1.17  -0.29 1.00      852      970
##    m1[35]   1.31   1.31 0.25 0.25   0.88   1.73 1.00      930     1139
##    m1[36]   0.92   0.92 0.25 0.24   0.50   1.33 1.00     1054     1346
##    m1[37]  -1.68  -1.67 0.28 0.28  -2.15  -1.23 1.00     2495     1731
##    m1[38]   0.63   0.63 0.24 0.23   0.24   1.01 1.00     1196     1277
##    m1[39]   4.06   4.06 0.31 0.31   3.54   4.56 1.00     1193     1533
##    m1[40]   3.76   3.76 0.33 0.32   3.22   4.31 1.00     1631     1711
##    m1[41]  -0.44  -0.43 0.26 0.25  -0.87  -0.03 1.00      844      953
##    m1[42]   2.48   2.48 0.28 0.27   2.03   2.94 1.00     1340     1623
##    m1[43]   3.24   3.23 0.30 0.30   2.73   3.75 1.00     1265     1324
##    m1[44]   1.39   1.39 0.23 0.21   1.02   1.75 1.00     1694     1386
##    m1[45]  -0.58  -0.57 0.24 0.24  -0.98  -0.18 1.00     2223     1639
##    m1[46]   2.34   2.34 0.24 0.24   1.94   2.74 1.00     1587     1350
##    m1[47]   3.38   3.38 0.29 0.28   2.91   3.84 1.00     1151     1482
##    m1[48]  -1.54  -1.53 0.28 0.27  -2.00  -1.10 1.00     2464     1790
##    m1[49]   0.00   0.00 0.23 0.23  -0.39   0.37 1.00     2050     1618
##    m1[50]   2.12   2.12 0.27 0.26   1.68   2.57 1.00     1047     1271

fn(f2)
## Initial log joint probability = -287.166 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       36      -18.3444   5.62153e-05    0.00205695      0.9495      0.9495       42    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -21.67 -21.38 1.81 1.69 -25.02 -19.37 1.01      804     1047
##    b[1]    -0.24  -0.25 0.27 0.26  -0.69   0.20 1.01      642     1004
##    b[2]     1.23   1.23 0.08 0.08   1.09   1.37 1.00     2428     1425
##    b[3]     1.43   1.45 0.37 0.36   0.83   2.01 1.01      527      951
##    b[4]     0.81   0.81 0.40 0.39   0.18   1.48 1.01      558      904
##    b[5]    -0.68  -0.68 0.53 0.54  -1.57   0.19 1.02      428      753
##    s        0.95   0.94 0.10 0.10   0.80   1.12 1.00     1989     1462
##    y1[1]   -0.69  -0.69 1.02 0.97  -2.41   0.94 1.00     1848     1954
##    y1[2]   -2.47  -2.48 1.05 1.03  -4.19  -0.78 1.00     1730     2029
##    y1[3]   -1.29  -1.28 0.98 0.96  -2.89   0.34 1.00     1875     1767
##    y1[4]   -0.31  -0.30 0.97 0.98  -1.94   1.25 1.00     2050     1929
##    y1[5]   -0.74  -0.72 1.00 0.98  -2.41   0.82 1.00     1938     2012
##    y1[6]    2.84   2.86 1.02 1.02   1.20   4.52 1.00     1910     1678
##    y1[7]   -1.60  -1.59 0.99 0.97  -3.26   0.06 1.00     2020     1888
##    y1[8]    2.97   3.00 0.99 0.97   1.32   4.54 1.00     1782     2003
##    y1[9]    2.22   2.20 0.99 0.98   0.56   3.86 1.00     1951     1820
##    y1[10]   0.81   0.78 0.98 0.96  -0.78   2.45 1.00     2198     2008
##    y1[11]  -0.83  -0.82 1.00 0.97  -2.49   0.78 1.00     2129     2148
##    y1[12]   1.82   1.81 1.00 0.97   0.22   3.47 1.00     1823     1810
##    y1[13]   4.52   4.54 1.02 1.05   2.82   6.21 1.00     1728     1804
##    y1[14]   1.86   1.88 0.97 0.93   0.31   3.44 1.00     1926     1995
##    y1[15]   1.09   1.05 0.99 0.97  -0.51   2.73 1.00     1844     1945
##    y1[16]   0.74   0.74 0.97 0.95  -0.84   2.39 1.00     2092     1905
##    y1[17]   3.55   3.57 1.00 0.97   1.90   5.23 1.00     1972     1918
##    y1[18]   3.12   3.12 1.00 1.00   1.48   4.77 1.00     1886     1815
##    y1[19]  -1.13  -1.11 1.03 0.99  -2.90   0.55 1.00     1994     1774
##    y1[20]  -2.15  -2.14 1.01 0.98  -3.83  -0.45 1.00     2118     1680
##    y1[21]  -3.71  -3.74 1.04 1.04  -5.36  -1.98 1.00     1601     1902
##    y1[22]   2.04   2.03 0.96 0.94   0.51   3.66 1.00     1905     1810
##    y1[23]   3.54   3.55 1.00 1.01   1.90   5.14 1.00     1943     1857
##    y1[24]   3.46   3.45 1.01 1.01   1.80   5.16 1.00     2011     1834
##    y1[25]   2.23   2.20 0.98 0.96   0.68   3.88 1.00     1925     2040
##    y1[26]   4.83   4.82 1.01 1.04   3.16   6.53 1.00     1968     1857
##    y1[27]  -1.39  -1.38 0.99 0.97  -3.05   0.24 1.00     2144     1878
##    y1[28]   0.94   0.95 1.01 0.99  -0.75   2.55 1.00     1945     1968
##    y1[29]  -1.50  -1.50 0.98 0.98  -3.09   0.12 1.00     1978     1934
##    y1[30]  -1.56  -1.56 1.02 1.00  -3.21   0.13 1.00     1847     1874
##    y1[31]   1.22   1.25 0.99 0.95  -0.46   2.77 1.00     1839     1911
##    y1[32]  -0.58  -0.56 1.00 1.01  -2.21   1.04 1.00     1895     1811
##    y1[33]  -0.59  -0.61 1.00 0.97  -2.23   1.03 1.00     1745     1883
##    y1[34]  -0.90  -0.90 0.97 0.94  -2.49   0.71 1.00     1987     1892
##    y1[35]   1.11   1.13 0.96 0.96  -0.50   2.64 1.00     1654     1929
##    y1[36]   1.13   1.15 0.97 0.96  -0.46   2.68 1.00     1918     2081
##    y1[37]  -1.79  -1.76 0.99 0.97  -3.45  -0.19 1.00     1917     1932
##    y1[38]   0.81   0.79 0.98 0.96  -0.82   2.41 1.00     1770     1767
##    y1[39]   4.20   4.19 0.99 0.98   2.58   5.81 1.00     2060     1675
##    y1[40]   3.95   3.94 1.03 1.02   2.27   5.69 1.00     1966     1629
##    y1[41]  -0.62  -0.63 0.99 0.97  -2.28   1.05 1.00     1886     1853
##    y1[42]   2.67   2.68 1.01 1.03   0.98   4.31 1.00     2018     1931
##    y1[43]   3.02   3.01 0.99 0.96   1.36   4.62 1.00     1945     2014
##    y1[44]   1.24   1.25 0.98 0.97  -0.39   2.84 1.00     1955     1906
##    y1[45]  -0.68  -0.67 0.94 0.93  -2.26   0.82 1.00     2112     1880
##    y1[46]   2.20   2.19 1.00 0.96   0.51   3.83 1.00     1894     1866
##    y1[47]   3.53   3.52 1.02 1.00   1.86   5.20 1.00     2037     1932
##    y1[48]  -1.66  -1.66 0.98 0.98  -3.24   0.00 1.00     1837     1730
##    y1[49]  -0.11  -0.09 0.99 0.97  -1.81   1.46 1.00     2107     2015
##    y1[50]   1.90   1.91 0.99 1.00   0.31   3.55 1.00     1966     1890
##    m1[1]   -0.68  -0.68 0.27 0.28  -1.14  -0.22 1.01      649     1003
##    m1[2]   -2.47  -2.47 0.38 0.37  -3.07  -1.84 1.00     1440     1542
##    m1[3]   -1.32  -1.32 0.33 0.33  -1.84  -0.78 1.00     1315     1329
##    m1[4]   -0.29  -0.29 0.31 0.31  -0.81   0.23 1.00     1176     1337
##    m1[5]   -0.73  -0.73 0.32 0.32  -1.26  -0.20 1.00     1208     1396
##    m1[6]    2.86   2.86 0.28 0.28   2.40   3.31 1.01     1209     1398
##    m1[7]   -1.62  -1.63 0.30 0.28  -2.11  -1.11 1.00     1896     1398
##    m1[8]    2.96   2.95 0.30 0.31   2.49   3.44 1.00     1580     1477
##    m1[9]    2.22   2.22 0.30 0.30   1.71   2.73 1.01     1172     1295
##    m1[10]   0.85   0.84 0.30 0.29   0.35   1.35 1.01     1149     1177
##    m1[11]  -0.83  -0.83 0.28 0.28  -1.29  -0.37 1.00      655     1002
##    m1[12]   1.84   1.84 0.29 0.28   1.36   2.32 1.00      745     1277
##    m1[13]   4.53   4.52 0.37 0.39   3.94   5.13 1.00     1683     1539
##    m1[14]   1.87   1.88 0.26 0.26   1.43   2.30 1.01     1140     1403
##    m1[15]   1.11   1.11 0.29 0.28   0.62   1.61 1.01     1144     1271
##    m1[16]   0.71   0.71 0.26 0.25   0.31   1.16 1.00     1568     1497
##    m1[17]   3.54   3.53 0.33 0.34   3.02   4.07 1.00     1617     1434
##    m1[18]   3.08   3.07 0.32 0.32   2.54   3.62 1.01     1241     1308
##    m1[19]  -1.10  -1.10 0.32 0.32  -1.61  -0.58 1.00     1290     1316
##    m1[20]  -2.12  -2.13 0.32 0.30  -2.64  -1.58 1.00     1995     1315
##    m1[21]  -3.67  -3.67 0.38 0.38  -4.30  -3.04 1.00      897     1466
##    m1[22]   2.05   2.05 0.27 0.27   1.62   2.48 1.01     1149     1417
##    m1[23]   3.53   3.54 0.30 0.31   3.05   4.02 1.01     1302     1380
##    m1[24]   3.45   3.45 0.30 0.31   2.96   3.93 1.01     1291     1346
##    m1[25]   2.21   2.21 0.30 0.29   1.71   2.71 1.00      786     1220
##    m1[26]   4.87   4.86 0.39 0.40   4.25   5.51 1.00     1703     1491
##    m1[27]  -1.38  -1.39 0.29 0.27  -1.85  -0.89 1.00     1852     1384
##    m1[28]   0.94   0.94 0.27 0.26   0.50   1.41 1.00      673     1058
##    m1[29]  -1.54  -1.55 0.29 0.27  -2.02  -1.04 1.00     1886     1406
##    m1[30]  -1.55  -1.55 0.35 0.35  -2.13  -0.98 1.00     1283     1512
##    m1[31]   1.25   1.25 0.26 0.27   0.81   1.68 1.01     1127     1352
##    m1[32]  -0.59  -0.60 0.27 0.27  -1.05  -0.14 1.01      647      972
##    m1[33]  -0.61  -0.61 0.30 0.29  -1.09  -0.11 1.01     1236     1304
##    m1[34]  -0.87  -0.88 0.28 0.28  -1.34  -0.42 1.00      655     1003
##    m1[35]   1.13   1.13 0.27 0.27   0.69   1.60 1.00      685     1167
##    m1[36]   1.11   1.11 0.29 0.28   0.62   1.61 1.01     1144     1271
##    m1[37]  -1.78  -1.79 0.30 0.28  -2.28  -1.26 1.00     1913     1419
##    m1[38]   0.80   0.81 0.27 0.26   0.36   1.24 1.01     1135     1439
##    m1[39]   4.20   4.20 0.32 0.33   3.67   4.72 1.00     1417     1537
##    m1[40]   3.93   3.93 0.35 0.35   3.35   4.51 1.01     1329     1340
##    m1[41]  -0.60  -0.61 0.27 0.27  -1.06  -0.15 1.01      647      972
##    m1[42]   2.66   2.66 0.31 0.31   2.14   3.17 1.01     1207     1265
##    m1[43]   3.04   3.03 0.33 0.32   2.50   3.59 1.00      897     1296
##    m1[44]   1.25   1.25 0.26 0.25   0.83   1.69 1.00     1519     1528
##    m1[45]  -0.69  -0.70 0.27 0.25  -1.13  -0.23 1.00     1762     1497
##    m1[46]   2.20   2.20 0.28 0.29   1.75   2.64 1.00     1537     1595
##    m1[47]   3.53   3.53 0.30 0.31   3.04   4.02 1.01     1301     1380
##    m1[48]  -1.65  -1.66 0.30 0.28  -2.14  -1.14 1.00     1900     1420
##    m1[49]  -0.13  -0.13 0.26 0.25  -0.54   0.32 1.00     1678     1547
##    m1[50]   1.93   1.93 0.29 0.28   1.45   2.42 1.00      754     1254

tb=tibble(xa=runif(n,-2,2),xb=runif(n,-2,2),
          ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
          y=rnorm(n,xa+xb-xa*xb+(ca=='b')*2+(cb=='b')*2-(ca=='b')*(cb=='b'),1))

grid.arrange(qplot(data=tb,xa,y,col=ca),
             qplot(data=tb,xa,y,col=cb),
             qplot(data=tb,xb,y,col=ca),
             qplot(data=tb,xb,y,col=cb))

fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)

mle=mdl$optimize(data=data)
mle

mcmc=goMCMC(mdl,data)

seeMCMC(mcmc,exc='m',ch=F)

y1=mcmc$draws('y1')
smy=summary(y1)

tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95

grid.arrange(
  qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')

tb=arrange(tb,y1)
grid.arrange(
  qplot(data=tb,1:n,y,col=ca)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  qplot(data=tb,1:n,y,col=cb)+
    geom_point(aes(x=1:n,y=y1),tb,col='black')+
    geom_line(aes(x=1:n,y=l5),tb,col='gray')+
    geom_line(aes(x=1:n,y=u5),tb,col='gray'),
  ncol=2
)
}


f0=formula(y~xa+xb+ca+cb)
f1=formula(y~xa+xb+ca*cb)
f2=formula(y~xa*xb+ca*cb)

fn(f0)
## Initial log joint probability = -1547.39 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       23      -45.4497    0.00023767   0.000932955           1           1       30    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -48.25 -47.89 1.85 1.67 -51.79 -45.94 1.01      642     1094
##    b[1]     0.34   0.35 0.41 0.40  -0.33   1.01 1.01      913      999
##    b[2]     1.09   1.09 0.20 0.19   0.78   1.41 1.00     2680     1372
##    b[3]     1.14   1.14 0.24 0.24   0.75   1.53 1.00     2296     1334
##    b[4]     0.77   0.77 0.49 0.47  -0.06   1.53 1.01      877     1010
##    b[5]     1.65   1.64 0.48 0.45   0.88   2.50 1.01      800      959
##    s        1.63   1.61 0.18 0.17   1.36   1.95 1.00     1928     1676
##    y1[1]    0.09   0.11 1.70 1.71  -2.64   2.88 1.00     1754     1763
##    y1[2]    1.14   1.12 1.71 1.73  -1.57   3.86 1.00     1855     1683
##    y1[3]    0.14   0.16 1.69 1.66  -2.66   2.85 1.00     1917     2104
##    y1[4]    5.95   5.93 1.69 1.63   3.22   8.72 1.00     2023     1573
##    y1[5]   -1.87  -1.88 1.66 1.64  -4.59   0.83 1.00     1836     1890
##    y1[6]    3.29   3.31 1.71 1.67   0.46   6.05 1.00     1886     1681
##    y1[7]    4.34   4.37 1.72 1.69   1.51   7.25 1.00     2013     1857
##    y1[8]    2.84   2.85 1.70 1.66   0.07   5.58 1.00     2126     2102
##    y1[9]   -0.70  -0.72 1.73 1.72  -3.58   2.17 1.00     1944     1742
##    y1[10]   0.51   0.56 1.70 1.71  -2.33   3.19 1.00     1851     1916
##    y1[11]   2.82   2.86 1.68 1.67   0.05   5.52 1.00     1889     1582
##    y1[12]   4.93   4.93 1.71 1.68   2.16   7.71 1.00     2157     1853
##    y1[13]   0.23   0.23 1.72 1.72  -2.56   3.11 1.00     1997     1974
##    y1[14]   0.81   0.84 1.74 1.71  -2.02   3.65 1.00     2004     2053
##    y1[15]  -2.39  -2.42 1.78 1.74  -5.25   0.53 1.00     1957     1873
##    y1[16]  -1.16  -1.13 1.72 1.74  -4.01   1.66 1.00     1954     1989
##    y1[17]   3.55   3.56 1.74 1.71   0.57   6.34 1.00     1989     1784
##    y1[18]   1.37   1.31 1.67 1.67  -1.30   4.23 1.00     2001     1937
##    y1[19]   2.66   2.62 1.76 1.70  -0.20   5.60 1.00     2107     1856
##    y1[20]   3.61   3.60 1.71 1.72   0.84   6.48 1.00     2033     1930
##    y1[21]  -0.40  -0.35 1.77 1.74  -3.30   2.47 1.00     1889     2004
##    y1[22]   3.12   3.14 1.71 1.69   0.28   5.94 1.00     1675     1875
##    y1[23]   3.39   3.40 1.67 1.64   0.55   6.19 1.00     1932     1890
##    y1[24]   2.90   2.86 1.77 1.74   0.08   5.78 1.00     1648     1818
##    y1[25]   1.15   1.13 1.78 1.77  -1.74   4.07 1.00     1854     1686
##    y1[26]  -1.56  -1.55 1.73 1.65  -4.43   1.31 1.00     1969     1895
##    y1[27]   3.97   3.98 1.72 1.70   1.20   6.78 1.00     1897     1906
##    y1[28]   1.41   1.40 1.73 1.69  -1.41   4.32 1.00     1913     1855
##    y1[29]   1.22   1.21 1.71 1.67  -1.57   4.13 1.00     1811     1962
##    y1[30]  -2.31  -2.36 1.76 1.74  -5.10   0.54 1.00     1725     1523
##    y1[31]   1.96   2.00 1.72 1.70  -0.85   4.71 1.00     1818     1843
##    y1[32]   0.20   0.19 1.73 1.71  -2.62   3.02 1.00     1918     1934
##    y1[33]   0.41   0.42 1.73 1.78  -2.37   3.24 1.00     2052     1973
##    y1[34]   2.16   2.15 1.73 1.70  -0.65   5.09 1.00     1839     1921
##    y1[35]  -1.24  -1.24 1.67 1.67  -4.02   1.52 1.00     2028     1844
##    y1[36]   1.48   1.47 1.76 1.70  -1.39   4.36 1.00     1908     2004
##    y1[37]   1.10   1.07 1.67 1.69  -1.61   3.83 1.00     1671     1792
##    y1[38]  -1.20  -1.20 1.72 1.66  -4.15   1.53 1.00     1623     1807
##    y1[39]   2.09   2.07 1.67 1.65  -0.72   4.81 1.00     1983     1884
##    y1[40]   3.06   3.05 1.75 1.72   0.21   5.89 1.00     2188     1961
##    y1[41]  -1.94  -1.92 1.74 1.68  -4.88   0.80 1.00     2201     1991
##    y1[42]   3.77   3.82 1.70 1.64   0.91   6.53 1.00     1904     1678
##    y1[43]   1.97   2.01 1.72 1.67  -0.85   4.79 1.00     1951     1970
##    y1[44]  -0.94  -0.91 1.67 1.69  -3.69   1.76 1.00     2202     1917
##    y1[45]  -0.05  -0.04 1.65 1.64  -2.75   2.66 1.00     1912     1856
##    y1[46]  -0.04  -0.04 1.75 1.69  -2.88   2.82 1.00     1916     1704
##    y1[47]   3.54   3.54 1.67 1.66   0.75   6.20 1.00     2159     1933
##    y1[48]  -2.94  -2.95 1.74 1.76  -5.82  -0.06 1.00     1613     1874
##    y1[49]  -1.24  -1.23 1.71 1.67  -4.00   1.56 1.00     2070     1917
##    y1[50]   0.82   0.82 1.67 1.57  -1.94   3.53 1.00     1790     1969
##    m1[1]    0.09   0.09 0.49 0.50  -0.75   0.86 1.00     1370     1308
##    m1[2]    1.17   1.20 0.38 0.39   0.52   1.76 1.00     1219     1314
##    m1[3]    0.13   0.15 0.38 0.38  -0.50   0.74 1.00     1264     1411
##    m1[4]    5.96   5.95 0.59 0.59   5.01   6.93 1.01     2195     1494
##    m1[5]   -1.82  -1.83 0.49 0.47  -2.63  -1.03 1.00     2057     1437
##    m1[6]    3.33   3.34 0.40 0.40   2.68   4.00 1.00     1651     1463
##    m1[7]    4.37   4.38 0.67 0.66   3.26   5.46 1.00     1661     1608
##    m1[8]    2.84   2.84 0.54 0.53   1.97   3.71 1.00     1342     1430
##    m1[9]   -0.67  -0.68 0.51 0.51  -1.49   0.17 1.00     1290     1372
##    m1[10]   0.51   0.51 0.49 0.49  -0.29   1.29 1.00     1508     1634
##    m1[11]   2.79   2.81 0.47 0.48   2.01   3.53 1.00     1309     1537
##    m1[12]   4.90   4.90 0.51 0.50   4.06   5.74 1.01     1795     1469
##    m1[13]   0.24   0.24 0.47 0.47  -0.53   0.99 1.00     2029     1595
##    m1[14]   0.90   0.89 0.67 0.65  -0.17   2.02 1.00     1899     1401
##    m1[15]  -2.41  -2.40 0.53 0.51  -3.31  -1.57 1.00     2032     1458
##    m1[16]  -1.17  -1.15 0.45 0.45  -1.90  -0.45 1.00     1216     1281
##    m1[17]   3.57   3.58 0.53 0.55   2.68   4.40 1.00     1525     1645
##    m1[18]   1.34   1.35 0.41 0.42   0.67   1.99 1.00     1729     1336
##    m1[19]   2.65   2.64 0.64 0.63   1.61   3.71 1.00     1691     1382
##    m1[20]   3.56   3.57 0.50 0.49   2.75   4.39 1.00     1418     1459
##    m1[21]  -0.33  -0.34 0.64 0.63  -1.39   0.71 1.00     1899     1524
##    m1[22]   3.06   3.05 0.52 0.50   2.23   3.92 1.01     1258     1109
##    m1[23]   3.39   3.40 0.41 0.41   2.72   4.06 1.00     1511     1348
##    m1[24]   2.91   2.91 0.58 0.57   1.96   3.85 1.00     1219     1108
##    m1[25]   1.17   1.15 0.64 0.63   0.11   2.21 1.00     1423     1027
##    m1[26]  -1.51  -1.50 0.46 0.44  -2.27  -0.77 1.00     1681     1551
##    m1[27]   3.96   3.97 0.53 0.54   3.13   4.86 1.00     2621     1646
##    m1[28]   1.43   1.46 0.40 0.40   0.77   2.06 1.00     1333     1361
##    m1[29]   1.25   1.25 0.61 0.60   0.22   2.22 1.01     1115     1169
##    m1[30]  -2.31  -2.30 0.54 0.54  -3.19  -1.45 1.00     1186     1402
##    m1[31]   1.97   1.96 0.61 0.62   1.00   2.98 1.01     1155     1289
##    m1[32]   0.19   0.19 0.52 0.51  -0.66   1.03 1.00     1533     1662
##    m1[33]   0.47   0.45 0.53 0.52  -0.39   1.35 1.01     1219     1290
##    m1[34]   2.13   2.15 0.53 0.54   1.22   2.99 1.01     1121     1421
##    m1[35]  -1.28  -1.27 0.64 0.62  -2.36  -0.19 1.00     1459     1416
##    m1[36]   1.43   1.42 0.58 0.59   0.49   2.36 1.01     1456     1343
##    m1[37]   1.12   1.12 0.43 0.44   0.42   1.80 1.00     1483     1402
##    m1[38]  -1.26  -1.26 0.49 0.48  -2.07  -0.48 1.00     2176     1363
##    m1[39]   2.03   2.06 0.42 0.42   1.33   2.68 1.00     1187     1365
##    m1[40]   3.08   3.08 0.51 0.52   2.23   3.89 1.00     1366     1302
##    m1[41]  -1.91  -1.91 0.57 0.57  -2.87  -0.95 1.00     1572     1539
##    m1[42]   3.75   3.74 0.51 0.52   2.95   4.61 1.00     2544     1519
##    m1[43]   1.95   1.94 0.58 0.59   1.00   2.94 1.00     1635     1284
##    m1[44]  -0.97  -0.96 0.43 0.41  -1.68  -0.28 1.00     1709     1433
##    m1[45]  -0.07  -0.07 0.41 0.40  -0.75   0.60 1.01      901      975
##    m1[46]  -0.06  -0.07 0.59 0.58  -1.01   0.93 1.00     1414     1381
##    m1[47]   3.56   3.57 0.54 0.56   2.71   4.46 1.00     2650     1653
##    m1[48]  -2.95  -2.95 0.58 0.58  -3.87  -2.05 1.00     1339     1536
##    m1[49]  -1.21  -1.20 0.44 0.42  -1.93  -0.50 1.00     1725     1543
##    m1[50]   0.84   0.86 0.39 0.39   0.19   1.45 1.01     1089     1400

fn(f1)
## Initial log joint probability = -100.856 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       36      -44.7365   6.08998e-05   0.000532295           1           1       46    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -48.08 -47.68 2.05 1.94 -51.83 -45.50 1.00      644      919
##    b[1]     0.11   0.11 0.51 0.50  -0.71   0.95 1.01      514      928
##    b[2]     1.08   1.08 0.19 0.18   0.76   1.40 1.00     2210     1318
##    b[3]     1.25   1.24 0.25 0.25   0.85   1.67 1.00     1227     1494
##    b[4]     1.20   1.18 0.67 0.64   0.06   2.30 1.01      498      723
##    b[5]     2.34   2.34 0.85 0.85   0.94   3.74 1.01      439      577
##    b[6]    -1.06  -1.08 1.04 0.97  -2.69   0.74 1.01      435      642
##    s        1.62   1.61 0.18 0.17   1.37   1.94 1.00     1791     1427
##    y1[1]   -0.03   0.01 1.67 1.69  -2.84   2.69 1.00     1792     1970
##    y1[2]    1.43   1.49 1.68 1.63  -1.40   4.13 1.00     1879     1890
##    y1[3]    0.35   0.34 1.65 1.63  -2.31   3.08 1.00     2066     1929
##    y1[4]    5.85   5.87 1.70 1.65   3.01   8.59 1.00     2030     1971
##    y1[5]   -1.70  -1.69 1.70 1.73  -4.61   1.03 1.00     1948     1972
##    y1[6]    3.12   3.14 1.67 1.67   0.35   5.78 1.00     2064     2012
##    y1[7]    4.26   4.27 1.81 1.75   1.23   7.21 1.00     1883     1866
##    y1[8]    2.77   2.77 1.72 1.68  -0.02   5.58 1.00     2046     1972
##    y1[9]   -0.58  -0.63 1.72 1.67  -3.38   2.26 1.00     1904     1950
##    y1[10]   0.26   0.27 1.69 1.69  -2.53   3.04 1.00     1921     1486
##    y1[11]   3.08   3.05 1.72 1.69   0.38   6.01 1.00     2180     1864
##    y1[12]   4.84   4.87 1.72 1.72   2.01   7.61 1.00     2176     1634
##    y1[13]  -0.09  -0.06 1.70 1.64  -2.93   2.68 1.00     1711     1807
##    y1[14]   0.81   0.83 1.77 1.74  -2.08   3.66 1.00     2223     1974
##    y1[15]  -2.34  -2.43 1.72 1.68  -5.11   0.49 1.00     2020     1974
##    y1[16]  -1.39  -1.43 1.72 1.72  -4.17   1.46 1.00     1633     1841
##    y1[17]   3.87   3.87 1.71 1.65   0.94   6.72 1.00     1893     1770
##    y1[18]   1.08   1.10 1.71 1.68  -1.67   3.90 1.00     2014     2015
##    y1[19]   2.94   2.95 1.72 1.74   0.19   5.82 1.00     1873     1973
##    y1[20]   3.53   3.53 1.71 1.69   0.77   6.32 1.00     1925     1931
##    y1[21]  -0.42  -0.40 1.77 1.77  -3.37   2.38 1.00     1899     1825
##    y1[22]   3.40   3.42 1.80 1.79   0.39   6.37 1.00     1564     1845
##    y1[23]   3.23   3.19 1.74 1.70   0.36   6.13 1.00     2135     2015
##    y1[24]   2.72   2.67 1.75 1.74  -0.08   5.61 1.00     1948     1884
##    y1[25]   1.10   1.12 1.73 1.75  -1.81   3.86 1.00     2053     1952
##    y1[26]  -1.43  -1.45 1.66 1.62  -4.29   1.30 1.00     1957     1915
##    y1[27]   3.66   3.67 1.73 1.71   0.82   6.52 1.00     2069     1767
##    y1[28]   1.66   1.69 1.72 1.64  -1.10   4.53 1.00     2075     1912
##    y1[29]   0.90   0.92 1.80 1.82  -2.10   3.87 1.00     1575     1744
##    y1[30]  -2.69  -2.68 1.82 1.82  -5.63   0.25 1.00     1607     1778
##    y1[31]   2.52   2.51 1.85 1.79  -0.45   5.54 1.00     1685     1845
##    y1[32]  -0.05  -0.08 1.77 1.75  -2.93   2.86 1.00     1968     1919
##    y1[33]   0.78   0.81 1.75 1.71  -2.09   3.66 1.00     1950     1930
##    y1[34]   2.30   2.32 1.74 1.67  -0.66   5.16 1.00     1937     1855
##    y1[35]  -1.02  -1.01 1.74 1.77  -3.88   1.85 1.00     1841     1880
##    y1[36]   1.66   1.68 1.72 1.68  -1.20   4.57 1.00     1821     1390
##    y1[37]   0.96   0.95 1.66 1.69  -1.75   3.65 1.00     2029     1899
##    y1[38]  -1.08  -1.12 1.73 1.71  -3.90   1.77 1.00     1849     1758
##    y1[39]   2.27   2.29 1.70 1.72  -0.51   5.11 1.00     1918     1650
##    y1[40]   2.96   3.02 1.76 1.71   0.02   5.84 1.00     1915     1933
##    y1[41]  -1.92  -1.93 1.71 1.68  -4.70   0.85 1.00     2231     2057
##    y1[42]   3.46   3.44 1.75 1.64   0.59   6.40 1.00     1661     1934
##    y1[43]   2.00   2.00 1.76 1.71  -0.84   4.80 1.00     2039     1738
##    y1[44]  -0.78  -0.78 1.69 1.66  -3.50   1.99 1.00     1941     1697
##    y1[45]  -0.34  -0.38 1.76 1.71  -3.26   2.60 1.00     1448     1773
##    y1[46]   0.21   0.29 1.79 1.83  -2.78   3.06 1.00     1675     1887
##    y1[47]   3.26   3.25 1.73 1.74   0.51   6.13 1.00     2209     1916
##    y1[48]  -3.37  -3.33 1.78 1.71  -6.31  -0.51 1.00     1566     1931
##    y1[49]  -1.11  -1.11 1.69 1.67  -3.80   1.67 1.00     2057     2059
##    y1[50]   1.07   1.10 1.71 1.70  -1.67   3.83 1.00     1892     1837
##    m1[1]   -0.05  -0.03 0.53 0.51  -0.92   0.79 1.00     1094     1063
##    m1[2]    1.39   1.40 0.43 0.41   0.68   2.11 1.00     1340     1244
##    m1[3]    0.31   0.31 0.41 0.40  -0.35   0.98 1.00     1569     1384
##    m1[4]    5.88   5.87 0.59 0.60   4.95   6.84 1.00     2809     1845
##    m1[5]   -1.69  -1.69 0.49 0.49  -2.49  -0.88 1.00     2168     1366
##    m1[6]    3.15   3.14 0.43 0.42   2.44   3.87 1.00     2123     1609
##    m1[7]    4.32   4.30 0.68 0.67   3.20   5.43 1.00     1740     1627
##    m1[8]    2.74   2.73 0.57 0.56   1.78   3.67 1.00     1322     1478
##    m1[9]   -0.62  -0.62 0.50 0.49  -1.43   0.19 1.00     1939     1596
##    m1[10]   0.29   0.27 0.52 0.51  -0.57   1.15 1.00     1756     1483
##    m1[11]   3.04   3.05 0.52 0.50   2.20   3.89 1.00     1405     1458
##    m1[12]   4.82   4.81 0.51 0.52   4.00   5.67 1.00     2573     1510
##    m1[13]  -0.07  -0.07 0.55 0.53  -1.00   0.81 1.00     1451     1544
##    m1[14]   0.87   0.88 0.68 0.67  -0.30   1.95 1.00     1809     1477
##    m1[15]  -2.34  -2.35 0.52 0.51  -3.18  -1.50 1.00     2373     1410
##    m1[16]  -1.40  -1.41 0.55 0.54  -2.27  -0.53 1.01      593     1016
##    m1[17]   3.87   3.89 0.59 0.57   2.92   4.86 1.00     1390     1411
##    m1[18]   1.07   1.07 0.47 0.47   0.27   1.83 1.00     1516     1381
##    m1[19]   2.97   2.97 0.71 0.68   1.78   4.13 1.00     1179     1262
##    m1[20]   3.48   3.49 0.49 0.48   2.69   4.32 1.00     2175     1584
##    m1[21]  -0.41  -0.40 0.66 0.64  -1.51   0.67 1.00     1620     1430
##    m1[22]   3.50   3.52 0.65 0.61   2.42   4.57 1.00      716      926
##    m1[23]   3.23   3.23 0.43 0.44   2.54   3.94 1.00     2305     1679
##    m1[24]   2.71   2.70 0.64 0.65   1.66   3.77 1.01      775     1206
##    m1[25]   1.07   1.07 0.63 0.59   0.03   2.09 1.00     1915     1669
##    m1[26]  -1.41  -1.41 0.46 0.46  -2.13  -0.66 1.00     2115     1475
##    m1[27]   3.69   3.70 0.59 0.57   2.72   4.68 1.00     1853     1521
##    m1[28]   1.68   1.69 0.46 0.44   0.94   2.44 1.00     1260     1400
##    m1[29]   0.92   0.93 0.73 0.70  -0.27   2.12 1.01      608     1039
##    m1[30]  -2.68  -2.69 0.69 0.67  -3.81  -1.58 1.01      507      897
##    m1[31]   2.52   2.52 0.78 0.78   1.24   3.78 1.01      628      935
##    m1[32]  -0.03  -0.04 0.55 0.54  -0.94   0.88 1.00     1770     1569
##    m1[33]   0.82   0.84 0.61 0.60  -0.22   1.82 1.00      845     1303
##    m1[34]   2.27   2.28 0.54 0.53   1.39   3.12 1.00     1725     1563
##    m1[35]  -0.98  -0.96 0.68 0.67  -2.12   0.13 1.00     1143     1335
##    m1[36]   1.74   1.75 0.64 0.64   0.67   2.78 1.00     1045     1283
##    m1[37]   0.89   0.89 0.47 0.47   0.12   1.67 1.00     1690     1527
##    m1[38]  -1.07  -1.06 0.51 0.50  -1.92  -0.24 1.00     1808     1516
##    m1[39]   2.26   2.27 0.46 0.45   1.51   3.01 1.00     1371     1414
##    m1[40]   2.99   3.00 0.51 0.49   2.17   3.84 1.00     2097     1661
##    m1[41]  -1.91  -1.92 0.56 0.54  -2.82  -0.99 1.00     2115     1560
##    m1[42]   3.48   3.49 0.57 0.54   2.56   4.41 1.00     1821     1597
##    m1[43]   1.91   1.91 0.60 0.58   0.91   2.85 1.00     1645     1521
##    m1[44]  -0.81  -0.80 0.45 0.44  -1.53  -0.08 1.00     1869     1454
##    m1[45]  -0.33  -0.33 0.53 0.52  -1.18   0.55 1.01      490      946
##    m1[46]   0.23   0.24 0.64 0.63  -0.86   1.28 1.00     1100     1341
##    m1[47]   3.27   3.29 0.61 0.57   2.29   4.26 1.00     1690     1625
##    m1[48]  -3.33  -3.33 0.73 0.72  -4.53  -2.15 1.01      544      910
##    m1[49]  -1.07  -1.07 0.45 0.44  -1.77  -0.35 1.00     1979     1471
##    m1[50]   1.01   1.01 0.41 0.40   0.34   1.68 1.00     1540     1448

fn(f2)
## Initial log joint probability = -423.203 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       39      -21.4026   6.40578e-05   0.000170812           1           1       46    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -25.71 -25.33 2.16 2.04 -29.85 -22.94 1.00      690     1089
##    b[1]     0.22   0.23 0.30 0.31  -0.27   0.69 1.00      650      903
##    b[2]     0.97   0.97 0.12 0.12   0.77   1.17 1.00     2759     1450
##    b[3]     1.06   1.06 0.16 0.16   0.80   1.33 1.01     1353     1269
##    b[4]     1.59   1.59 0.40 0.40   0.93   2.24 1.00      562      713
##    b[5]     1.81   1.79 0.55 0.54   0.94   2.73 1.00      421      605
##    b[6]    -0.96  -0.96 0.12 0.12  -1.16  -0.76 1.00     2380     1685
##    b[7]    -1.04  -1.04 0.66 0.63  -2.12   0.03 1.00      384      571
##    s        1.03   1.02 0.12 0.11   0.87   1.24 1.00     2002     1679
##    y1[1]    0.81   0.80 1.11 1.10  -1.00   2.64 1.00     2046     1918
##    y1[2]    1.96   1.99 1.08 1.08   0.21   3.76 1.00     2025     1893
##    y1[3]    0.77   0.74 1.04 1.01  -0.95   2.47 1.00     1918     1998
##    y1[4]    3.65   3.64 1.14 1.13   1.76   5.52 1.00     1793     1860
##    y1[5]   -2.29  -2.29 1.10 1.06  -4.11  -0.46 1.00     2053     1860
##    y1[6]    3.18   3.20 1.06 1.08   1.35   4.86 1.00     1945     1971
##    y1[7]    0.77   0.76 1.23 1.23  -1.23   2.84 1.00     2126     1791
##    y1[8]    1.37   1.36 1.12 1.13  -0.44   3.20 1.00     1957     2013
##    y1[9]    0.09   0.09 1.10 1.13  -1.73   1.83 1.00     1968     1881
##    y1[10]  -0.31  -0.32 1.09 1.08  -2.09   1.48 1.00     2098     2012
##    y1[11]   2.78   2.78 1.14 1.08   0.92   4.63 1.00     1931     1863
##    y1[12]   3.65   3.63 1.07 1.06   1.93   5.40 1.00     1910     1895
##    y1[13]  -0.84  -0.81 1.09 1.03  -2.67   0.93 1.00     1695     1918
##    y1[14]   3.38   3.38 1.18 1.15   1.36   5.34 1.00     1939     1893
##    y1[15]  -3.78  -3.75 1.11 1.10  -5.59  -2.00 1.00     2259     2020
##    y1[16]  -1.29  -1.33 1.09 1.05  -3.07   0.56 1.00     1871     1720
##    y1[17]   2.90   2.90 1.06 1.07   1.23   4.65 1.00     1803     1751
##    y1[18]   1.04   1.05 1.09 1.10  -0.72   2.78 1.00     1666     1921
##    y1[19]   4.91   4.89 1.16 1.12   3.01   6.80 1.00     1933     1929
##    y1[20]   3.42   3.39 1.08 1.05   1.63   5.23 1.00     1965     1933
##    y1[21]   2.17   2.20 1.16 1.14   0.22   4.01 1.00     2070     1802
##    y1[22]   3.13   3.10 1.14 1.13   1.25   4.97 1.00     1739     1971
##    y1[23]   3.13   3.13 1.03 1.00   1.38   4.82 1.00     2011     2014
##    y1[24]   1.72   1.68 1.10 1.06  -0.07   3.56 1.00     1965     1799
##    y1[25]   2.12   2.08 1.10 1.07   0.38   3.96 1.00     2056     2040
##    y1[26]  -1.92  -1.95 1.04 1.03  -3.60  -0.18 1.00     1907     2015
##    y1[27]   5.19   5.22 1.13 1.10   3.31   7.07 1.00     1921     2145
##    y1[28]   2.20   2.20 1.10 1.06   0.37   3.97 1.00     1957     1869
##    y1[29]   2.31   2.34 1.14 1.11   0.39   4.16 1.00     1833     1775
##    y1[30]  -3.49  -3.44 1.09 1.06  -5.29  -1.70 1.00     1630     1729
##    y1[31]   2.64   2.65 1.17 1.24   0.76   4.48 1.00     1396     1818
##    y1[32]  -0.80  -0.80 1.14 1.09  -2.72   1.07 1.00     1917     1765
##    y1[33]   0.38   0.41 1.09 1.08  -1.41   2.14 1.00     1699     1891
##    y1[34]   3.39   3.41 1.09 1.05   1.56   5.19 1.00     2024     1910
##    y1[35]  -2.88  -2.86 1.18 1.14  -4.91  -0.93 1.00     1922     1921
##    y1[36]   2.76   2.76 1.11 1.07   0.98   4.58 1.00     1750     1919
##    y1[37]   0.60   0.58 1.05 1.01  -1.11   2.31 1.00     2019     2032
##    y1[38]  -0.73  -0.72 1.06 1.03  -2.50   1.00 1.00     2079     1875
##    y1[39]   2.52   2.50 1.10 1.06   0.79   4.34 1.00     1962     1884
##    y1[40]   3.32   3.32 1.08 1.10   1.56   5.05 1.00     1861     2034
##    y1[41]  -2.30  -2.33 1.10 1.08  -4.10  -0.47 1.00     1894     1893
##    y1[42]   4.78   4.80 1.08 1.07   2.98   6.56 1.00     2186     1921
##    y1[43]   2.33   2.33 1.09 1.05   0.54   4.15 1.00     1957     1973
##    y1[44]  -0.78  -0.75 1.10 1.09  -2.67   1.02 1.00     1779     1766
##    y1[45]  -0.20  -0.19 1.09 1.04  -2.03   1.57 1.00     2060     2010
##    y1[46]   0.02  -0.02 1.12 1.10  -1.77   1.93 1.00     1961     1800
##    y1[47]   5.15   5.15 1.12 1.10   3.33   6.96 1.00     1777     1856
##    y1[48]  -4.87  -4.88 1.12 1.12  -6.71  -3.02 1.00     1646     1777
##    y1[49]  -1.21  -1.25 1.09 1.07  -3.02   0.68 1.00     2111     2109
##    y1[50]   1.56   1.57 1.08 1.03  -0.23   3.27 1.00     1998     1890
##    m1[1]    0.83   0.83 0.35 0.33   0.25   1.39 1.00     1167     1178
##    m1[2]    1.90   1.90 0.28 0.28   1.44   2.35 1.00     1107     1263
##    m1[3]    0.76   0.76 0.26 0.26   0.33   1.20 1.00     1335     1310
##    m1[4]    3.62   3.61 0.47 0.47   2.85   4.41 1.00     2459     1349
##    m1[5]   -2.26  -2.27 0.31 0.31  -2.76  -1.76 1.00     2097     1527
##    m1[6]    3.17   3.18 0.27 0.27   2.74   3.62 1.00     1683     1735
##    m1[7]    0.77   0.76 0.61 0.61  -0.21   1.75 1.00     2135     1509
##    m1[8]    1.35   1.35 0.38 0.37   0.73   1.99 1.00     1495     1559
##    m1[9]    0.10   0.10 0.32 0.32  -0.43   0.62 1.00     1673     1472
##    m1[10]  -0.30  -0.30 0.34 0.32  -0.87   0.25 1.00     1704     1643
##    m1[11]   2.81   2.81 0.33 0.34   2.25   3.34 1.00     1282     1311
##    m1[12]   3.65   3.66 0.36 0.35   3.08   4.25 1.00     2231     1573
##    m1[13]  -0.83  -0.83 0.36 0.35  -1.43  -0.27 1.00     1426     1481
##    m1[14]   3.35   3.36 0.54 0.52   2.48   4.21 1.00     2214     1608
##    m1[15]  -3.75  -3.75 0.37 0.37  -4.37  -3.15 1.00     2540     1639
##    m1[16]  -1.30  -1.30 0.33 0.32  -1.85  -0.77 1.00      741      943
##    m1[17]   2.88   2.88 0.40 0.39   2.24   3.52 1.00     1359     1430
##    m1[18]   1.04   1.04 0.29 0.28   0.56   1.51 1.00     1324     1445
##    m1[19]   4.91   4.90 0.50 0.49   4.09   5.70 1.00     1666     1654
##    m1[20]   3.43   3.43 0.31 0.31   2.93   3.96 1.00     2078     1713
##    m1[21]   2.15   2.16 0.54 0.52   1.29   3.02 1.00     1989     1588
##    m1[22]   3.13   3.13 0.43 0.42   2.46   3.84 1.00      762      889
##    m1[23]   3.10   3.10 0.27 0.27   2.66   3.54 1.00     1871     1815
##    m1[24]   1.71   1.71 0.40 0.40   1.04   2.35 1.00     1098     1362
##    m1[25]   2.15   2.14 0.43 0.40   1.47   2.87 1.00     2401     1589
##    m1[26]  -1.88  -1.88 0.29 0.29  -2.35  -1.41 1.00     2023     1409
##    m1[27]   5.13   5.14 0.39 0.39   4.51   5.76 1.00     1703     1474
##    m1[28]   2.21   2.21 0.30 0.31   1.72   2.69 1.00      985     1213
##    m1[29]   2.35   2.35 0.47 0.47   1.57   3.09 1.00      813     1196
##    m1[30]  -3.47  -3.46 0.42 0.42  -4.17  -2.78 1.00      675     1020
##    m1[31]   2.64   2.63 0.51 0.52   1.84   3.48 1.00      639      909
##    m1[32]  -0.79  -0.79 0.36 0.35  -1.39  -0.20 1.00     1792     1718
##    m1[33]   0.38   0.38 0.40 0.40  -0.27   1.06 1.00      910     1076
##    m1[34]   3.37   3.37 0.37 0.35   2.77   3.97 1.00     1586     1421
##    m1[35]  -2.89  -2.89 0.51 0.49  -3.73  -2.05 1.00     1254     1347
##    m1[36]   2.78   2.76 0.42 0.43   2.10   3.44 1.00     1318     1653
##    m1[37]   0.59   0.60 0.30 0.28   0.08   1.08 1.00     1532     1712
##    m1[38]  -0.74  -0.73 0.32 0.31  -1.26  -0.22 1.00     1420     1725
##    m1[39]   2.48   2.48 0.30 0.30   1.99   2.96 1.00     1224     1256
##    m1[40]   3.29   3.29 0.33 0.32   2.77   3.84 1.00     2119     1452
##    m1[41]  -2.27  -2.27 0.35 0.36  -2.84  -1.69 1.00     2052     1722
##    m1[42]   4.77   4.78 0.37 0.37   4.17   5.37 1.00     1657     1451
##    m1[43]   2.36   2.36 0.38 0.37   1.74   2.97 1.01     1737     1558
##    m1[44]  -0.73  -0.73 0.28 0.27  -1.18  -0.27 1.00     1594     1644
##    m1[45]  -0.19  -0.18 0.31 0.32  -0.70   0.30 1.00      616      792
##    m1[46]  -0.01   0.01 0.41 0.41  -0.69   0.66 1.00     1183     1323
##    m1[47]   5.14   5.15 0.43 0.43   4.46   5.83 1.00     1639     1394
##    m1[48]  -4.85  -4.85 0.47 0.47  -5.63  -4.07 1.00      793     1255
##    m1[49]  -1.22  -1.21 0.28 0.28  -1.67  -0.75 1.00     1793     1489
##    m1[50]   1.57   1.58 0.27 0.27   1.14   2.02 1.00     1328     1070